Code to load tools and data
 # packages
    require(arm)
    require(data.table)
    require(effects)
    require(foreach)
    require(ggimage)
    require(ggplot2)
    require(ggpubr)
    require(ggsci)
    require(ggtext)
    require(grid)
    require(gtable)
    require(here)
    require(kableExtra)
    require(MASS)
    require(multcomp)
    require(optimx)
    require(performance)  
    require(PerformanceAnalytics)
    require(png)
    require(RColorBrewer)
    require(rmeta)
    require(rphylopic)
    require(scales)
    require(viridis)
 # constants
    save_plot = TRUE
    round_ = 3 # number of decimal places to round model coefficients
    nsim = 5000 # number of simulations to extract estimates and 95%CrI
    ax_lines = "grey60" # defines color of the axis lines
    #colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
    set.seed(42)
    #width_ = .7 # spacing between error bars
    #col_ = c(brewer.pal(n =12, name = "Paired"), 'grey30','grey80')
 # functions
    # to add images to panels
      annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data) {
        layer(data = data, stat = StatIdentity, position = PositionIdentity, 
            geom = ggplot2:::GeomCustomAnn,
            inherit.aes = TRUE, params = list(grob = grob, 
                                              xmin = xmin, xmax = xmax, 
                                              ymin = ymin, ymax = ymax))
      }
    # to remove ggplot components
        gtable_filter_remove <- function (x, name, trim = TRUE){
          matches <- !(x$layout$name %in% name)
          x$layout <- x$layout[matches, , drop = FALSE]
          x$grobs <- x$grobs[matches]
          if (trim) 
            x <- gtable_trim(x)
          x
        }
    # customized ggplot theme
        theme_MB = theme(  
                  title = element_text(size=8, colour="grey30"),
                  axis.line = element_blank(),
                  #axis.line = element_line(colour="grey70", size=0.25),
                  axis.title = element_text(size=7, colour="grey30"),
                  axis.title.y = element_text(vjust=3.5),
                  axis.title.x = element_text(vjust=1),
                  axis.text = element_text(size=6),#, vjust = 0.5, hjust=1),# margin=units(0.5,"mm")),
                  axis.ticks.length=unit(0.5,"mm"),
                  axis.ticks = element_line(colour = "grey70", size = 0.1),
                  #axis.ticks.margin,
                  
                  strip.text.x = element_text(size = 6, color="grey30",  margin=margin(1,1,1,1,"mm")), #grey50
                  strip.text.y = element_text(size = 6, color="grey30",  margin=margin(1,1,1,1,"mm")), #grey50
                  strip.background = element_rect(fill="grey99",colour="grey70", size=0.25),
                    #strip.background = element_blank(), 
                    #strip.text = element_blank(),
                  panel.spacing = unit(0, "mm"),
                  panel.background=element_blank(),
                  panel.border = element_rect(colour="grey70", size=0.1, fill = NA), #panel.border=element_blank(),
                  panel.grid = element_blank(),

                  legend.text=element_text(size=6),
                  legend.title=element_text(size=6),
                  legend.key = element_rect(colour = NA, fill = NA),
                  legend.key.height= unit(0.5,"line"),
                  legend.key.width = unit(0.25, "cm"),
                  legend.margin = margin(0,0,0,0, unit="cm"),
                  legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
                  legend.background = element_blank()
                  )  
    # for estimates
        est_out =function(model = m, label = "", nsim = 5000){
            bsim = sim(model, n.sim=5000) 
            v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
            sd = apply(bsim@fixef, 2, sd)
            data.table(predictor=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,], sd = sd, model = paste(label, "N =", nobs(model)))
          }
    # change color
      change_col = function(replace_black, theimg) {
        r_b = col2rgb(replace_black) / 255
        #theimg[theimg == 1] <- 2
        for (i in 1:3) {
            theimg[,,i][theimg[,,i] == 0] <- r_b[i]
        }
        return(theimg)
        }
    # for Supplementary Table output based on sim
      m_out = function(model = m, type = "mixed", 
        name = "define", dep = "define", fam = 'Gaussian',
        round_ = 3, nsim = 5000, aic = FALSE, save_sim = here::here('Data/model_sim/'), back_tran = FALSE, perc_ = 1){
          # perc_ 1 = proportion or 100%
        bsim = sim(model, n.sim=nsim)  
        
        if(save_sim!=FALSE){save(bsim, file = paste0(save_sim, name,'.RData'))}
       
        if(type != "mixed"){
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 

          if(back_tran == TRUE & fam == "binomial"){
           v = perc_*plogis(v)
           ci = perc_*plogis(ci)
           }
          if(back_tran == TRUE & fam == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam == "Poisson"){
           v = exp(v)
           ci = exp(ci)
          }

         oi=data.frame(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
          rownames(oi) = NULL
          oi$estimate_r=round(oi$estimate,round_)
          oi$lwr_r=round(oi$lwr,round_)
          oi$upr_r=round(oi$upr,round_)
          if(perc_ == 100){
           oi$estimate_r = paste0(oi$estimate_r,"%")
           oi$lwr_r = paste0(oi$lwr_r,"%")
           oi$upr_r = paste0(oi$upr_r,"%")
          }
         x=data.table(oi[c('type',"effect", "estimate_r","lwr_r",'upr_r')]) 
       
        }else{
         v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
         ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 

         if(back_tran == TRUE & fam == "binomial"){
          v = perc_*plogis(v)
          ci = perc_*plogis(ci)
         }
          if(back_tran == TRUE & fam == "binomial_logExp"){
                v = perc_*(1-plogis(v))
                ci = perc_*(1-plogis(ci))
                ci = rbind(ci[2,],ci[1,])
               }

          if(back_tran == TRUE & fam == "Poisson"){
            v = exp(v)
            ci = exp(ci)
         }

        oi=data.table(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
            rownames(oi) = NULL
            oi[,estimate_r := round(estimate,round_)]
            oi[,lwr_r := round(lwr,round_)]
            oi[,upr_r :=round(upr,round_)]
            if(perc_ == 100){
             oi[,estimate_r := paste0(estimate_r,"%")]
             oi[,lwr_r := paste0(lwr_r,"%")]
             oi[,upr_r := paste0(upr_r,"%")]
            }
         oii=oi[,c('type',"effect", "estimate_r","lwr_r",'upr_r')] 
        
         l=data.frame(summary(model)$varcor)
         l = l[is.na(l$var2),]
         l$var1 = ifelse(is.na(l$var1),"",l$var1)
         l$pred = paste(l$grp,l$var1)

         q050={}
         q025={}
         q975={}
         pred={}
         
         # variance of random effects
         for (ran in names(bsim@ranef)) {
           ran_type = l$var1[l$grp == ran]
           for(i in ran_type){
            q050=c(q050,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.5)))
            q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.025)))
            q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.975)))
            pred= c(pred,paste(ran, i))
            }
           }
         # residual variance
         q050=c(q050,quantile(bsim@sigma^2, prob=c(0.5)))
         q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
         q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
         pred= c(pred,'Residual')

         ri=data.table(type='random',effect=pred, estimate_r=round(100*q050/sum(q050)), lwr_r=round(100*q025/sum(q025)), upr_r=round(100*q975/sum(q975)))
           
         ri[lwr_r>upr_r, lwr_rt := upr_r]
         ri[lwr_r>upr_r, upr_rt := lwr_r]
         ri[!is.na(lwr_rt), lwr_r := lwr_rt]
         ri[!is.na(upr_rt), upr_r := upr_rt]
         ri$lwr_rt = ri$upr_rt = NULL

         ri[,estimate_r := paste0(estimate_r,'%')]
         ri[,lwr_r := paste0(lwr_r,'%')]
         ri[,upr_r := paste0(upr_r,'%')]
        
        x = data.table(rbind(oii,ri))
        }
        
        x[1, model := name]                                                                
        x[1, response := dep]                                                                
        x[1, error_structure := fam]      
        N = length(resid(model))                                                          
        x[1, N := N ]                                                                

        x=x[ , c('model', 'response', 'error_structure', 'N', 'type',"effect", "estimate_r","lwr_r",'upr_r')] 

        if (aic == TRUE){   
            x[1, AIC := AIC(update(model,REML = FALSE))] 
            }
        if (aic == "AICc"){
            aicc = AICc(model)
            x[1, AICc := aicc] 
        }
        if(type == "mixed" & nrow(x[type=='random' & estimate_r =='0%'])==0){
          R2_mar = as.numeric(r2_nakagawa(model)$R2_marginal)
          R2_con = as.numeric(r2_nakagawa(model)$R2_conditional)
          x[1, R2_mar := R2_mar]
          x[1, R2_con := R2_con]
         }
        x[is.na(x)] = ""
        return(x)
      } 
    # model assumption function
      m_ass = function(name = 'define', mo = m0, dat = d, fixed = NULL, categ = NULL, trans = "none", spatial = TRUE, temporal = TRUE, PNG = TRUE, outdir = 'outdir', n_col=8, width_ = 10, height_ = 5.5){
       l=data.frame(summary(mo)$varcor)
       l = l[is.na(l$var2),]
       nt = if(temporal==TRUE){1}else{0}
       ns = if(spatial==TRUE){7}else{0}
       n = 3+nrow(l)-1+length(fixed)+length(categ) +  nt +  ns
     
       if(PNG == TRUE){
        png(paste(outdir,name, ".png", sep=""), width=width_,height=height_,units="in",res=600) # width = 6
        par(mfrow=c(4, n_col),tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
            oma = c(1,1,2,1),
            mar = c(2, 2, 2, 1), mgp=c(1,0,0)
            )
         }else{
          dev.new(width=width_,height=height_)
          par(mfrow=c(4,n_col), tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
            oma = c(1,1,2,1),
            mar = c(2, 2, 2, 1), mgp=c(1,0,0)
            )
        }
       
       scatter.smooth(fitted(mo),resid(mo),col='grey');abline(h=0, lty=2, col ='red')
       scatter.smooth(fitted(mo),sqrt(abs(resid(mo))), col='grey')
       qqnorm(resid(mo), main=list("Normal Q-Q Plot: residuals"),col='grey');qqline(resid(mo), col = 'red')
       #unique(l$grp[l$grp!="Residual"])
       for(i in unique(l$grp[l$grp!="Residual"])){
        #i = "mean_year"
        ll=ranef(mo)[names(ranef(mo))==i][[1]]
        if(ncol(ll)==1){
         qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey',);qqline(ll[,1], col ='red')
         }else{
          qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey');qqline(ll[,1], col ='red')
          qqnorm(ll[,2], main = paste(i,names(ll)[2]),col='grey');qqline(ll[,2], col ='red')
         }
        }
        
       # variables
         scatter={} 
         for (i in rownames(summary(mo)$coef)) {
              #i = "lat_abs"
            j=sub("\\).*", "", sub(".*\\(", "",i)) 
            scatter[length(scatter)+1]=j
          }
          x = data.frame(scatter=unique(scatter)[2:length(unique(scatter))],
                          log_ = grepl("log",rownames(summary(mo)$coef)[2:length(unique(scatter))]), stringsAsFactors = FALSE)
          for (i in 1:length(fixed)){
              jj =fixed[i]
              variable=dat[, ..jj][[1]]
              if(trans[i]=='log'){
              scatter.smooth(resid(mo)~log(variable),xlab=paste('log(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='abs'){
              scatter.smooth(resid(mo)~abs(variable),xlab=paste('abs(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='sin'){scatter.smooth(resid(mo)~sin(variable),xlab=paste('sin(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else if(trans[i]=='cos'){scatter.smooth(resid(mo)~cos(variable),xlab=paste('cos(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
              }else{
              scatter.smooth(resid(mo)~variable,xlab=jj,col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
            }
           }
          
          if(length(categ)>0){
            for(i in categ){
               variable=dat[, ..i][[1]]
                boxplot(resid(mo)~variable, medcol='grey', whiskcol='grey', staplecol='grey', boxcol='grey', outcol='grey', xlab = i);abline(h=0, lty=3, lwd=1, col = 'red')
               }
          }     
              
        if(temporal == TRUE){
            acf(resid(mo), type="p", main=list("Temporal autocorrelation:\npartial series residual"))
            }
        if(spatial == TRUE){    
          spdata=data.frame(resid=resid(mo), x=dat$Lon, y=dat$Lat)
            spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x, spdata$y,col=spdata$col, cex=as.numeric(spdata$cex), pch= 16, main=list('Spatial distribution of residuals', cex=0.8), xlab = 'Longitude', ylab = 'Latitude')
            legend("topright", pch=16, legend=c('>0','<0'), ,col=c(rgb(83,95,124,100, maxColorValue = 255),rgb(253,184,19,100, maxColorValue = 255)))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('residual >=0'), xlab = 'Longitude', ylab = 'Latitude')

          # EU
          dat$res = resid(mo)
          spdata=data.frame(resid = dat$res[dat$Country!='Australia'], x=dat$Lon[dat$Country!='Australia'], y=dat$Lat[dat$Country!='Australia'])
          spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('EU -  residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('EU residuals >=0)'), xlab = 'Longitude', ylab = 'Latitude')

          # Australia
          spdata=data.frame(resid = dat$res[dat$Country=='Australia'], x=dat$Lon[dat$Country=='Australia'], y=dat$Lat[dat$Country=='Australia'])
          spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
            #cex_=c(1,2,3,3.5,4)
            cex_=c(1,1.5,2,2.5,3)
            spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
          plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('Australia residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
          plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('Australia residuals >=0'), xlab = 'Longitude', ylab = 'Latitude')
          }
       
       mtext(stringr::str_wrap(paste(paste0(name," model: "), slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''), width = ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)+10), side = 3, line = 0, cex=0.5,outer = TRUE, col = 'darkblue') #ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)
       if(PNG==TRUE){dev.off()}
      }
    
 # data
    ph  =  fread(here::here('Data/phylopic.txt'))
    setnames(ph, old = c('Name', 'Code'), new = c('genus2', 'uid'))

    t = fread(here::here('Data/taxonomy.txt'))
    
    g = fread(here::here('Data/google_mobility.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
    g[, Year := as.integer(substring(date, nchar(date)-3, nchar(date)))]
    g[nchar(date)==9, date:=paste0('0',date)]
    g[, date_ :=as.Date(date, format = '%d.%m.%Y')]
    g[, Day :=yday(date_)]
    g[country_region!='Australia', Day := Day-92 +1] # 1 April = start of breeding season (1st day) = 92 day of the year 
    g[country_region=='Australia', Day := Day-228 +1] # 15 Augusst = start of breeding season (1st day) = 228 day of the year 
    setnames(g, old = 'country_region', new ='Country')

    d = fread(here::here('Data/data_corrected.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
    # add data and weekdays
    x = fread(here::here("Data/date.txt"))[, .(IDObs, Date_corr)]
    x[, date_:=as.Date(Date_corr, "%d.%m.%Y" )]
    x[, weekday := weekdays(date_)]
    d =  merge(d, x[, .(IDObs, date_, weekday)], by = "IDObs")
    #d[, date_ := as.Date(Day, origin = paste0(Year, '-01-01'))]

    # adjust correct assignment of season (Year) for Australia
    d[Country == 'Australia' & Year == 2020 & Covid == 0, Year:=2019]
    d[Country == 'Australia' & Year == 2021 & Day>139, Year:=2020]
    d[Country == 'Australia' & Year == 2022 & Day>139, Year:=2021]

    d = d[order(Year, IDLocality, Day, Hour)]
    d[Country %in% c("Czech_Republic", "Czech Republic"), Country := "Czechia"]
    d[, genus := sub("_.*", "", Species)]
    d[, sp_day_year := paste(Year, Species, Day, sep="_")]
    d[, sp_loc := paste(Species, IDLocality, sep="_")]
    d[, sp_country := paste(Species, Country, sep="\n")]
    d[, rad:=(2*pi*Hour) / 24]
    d[, Day_:= Day]
    #d[, FID_z := scale(FID), by = Species]
    #d[, SD_z := scale(SD), by = Species]
    d[, FID_ln := log(FID)]
    d[, SD_ln := log(SD)]
    d[, body_ln := log(BodyMass)]
    d[, flock_ln := log(FlockSize)]
    d[, weekday := weekdays(date_)]
    #d[Country == 'Australia', Day_:= abs(Day - 189)]
    
    d1 = d[Covid == 1, .N, by = Species]
    d2 = d[Covid == 0, .N, by = Species]
    setnames(d1, old = 'N', new ='N_during')
    setnames(d2, old = 'N', new ='N_before')
    dd = merge(d1,d2) # species with data before and during
    da = merge(d1,d2, all = TRUE)

    d1p = d[Country!='Poland' & Covid == 1, .N, by = Species]
    d2p = d[Country!='Poland' & Covid == 0, .N, by = Species]
    setnames(d1p, old = 'N', new ='N_during')
    setnames(d2p, old = 'N', new ='N_before')
    ddp = merge(d1p,d2p) # species with data before and during, but without Poland data

    p1 = d[Covid == 1, .N, by = .(IDLocality, Species)]
    p2 = d[Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1, old = 'N', new ='N_during')
    setnames(p2, old = 'N', new ='N_before')
    pp = merge(p1,p2)  # species-localities with data before and during
    pa = merge(p1,p2, all = TRUE)

    p1p = d[Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
    p2p = d[Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1p, old = 'N', new ='N_during')
    setnames(p2p, old = 'N', new ='N_before')
    ppp = merge(p1p,p2p)  # species-localities with data before and during,but without Poland data 
   
    p1p4 = d[Year!=2014 & Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
    p2p4 = d[Year!=2014 &Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
    setnames(p1p4, old = 'N', new ='N_during')
    setnames(p2p4, old = 'N', new ='N_before')
    ppp4 = merge(p1p4,p2p4) # species-localities with data before and during, but without Poland & 2014 data

    s = d[Covid == 1]
    s[, Nsp := .N, by ='Species']
    s[, sp := gsub('[_]', ' ', Species)]
    # add google mobility
    s = merge(s, g[,.(Country,  date_, parks_percent_change_from_baseline)], all.x = TRUE)
    s[, year_ := as.character(Year)]  

    ss = s[!is.na(parks_percent_change_from_baseline)]
    ss[, country_year := paste(Country, Year)] #table(paste(s$Country, s$Year))   
    ss[parks_percent_change_from_baseline<0, google := 'before_zero']
    ss[parks_percent_change_from_baseline>0, google := 'after_zero']
    ss[, sp_country_google:= paste(sp_country, google)]

    g[, weekday := weekdays(date_)]

Introduction

To facilitate transparency, the following document contains example code used to generate the results of the manuscript. Thus, apart from the supplementary figures and tables, below we display also main text figures. The figures and tables are order according to their appearance in the main text. The code is displayed upon clicking the “code” icon at the top right above each display item.


Effect sizes for Period

# predictions
# full
ms <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1| genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
data = d, REML = FALSE,
control <- lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_ms <- est_out(ms, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_ms[, control_for_starting_distance := "yes"]

mx <- lmer(scale(log(FID)) ~
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
    data = d, REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
        )  
# (Covid|IDLocality) +
est_mx <- est_out(mx, "All: (1|Year) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
est_mx[, control_for_starting_distance := "no"]

# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  cs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = d[Country == "Czechia"], REML = FALSE
    )
  est_cs <- est_out(cs, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
  est_cs[, control_for_starting_distance := "yes"]
  
  cx <- lmer(scale(log(FID)) ~
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = d[Country == "Czechia"], REML = FALSE
    )
  est_cx <- est_out(cx, "Czechia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
  est_cx[, control_for_starting_distance := "no"]

# FI
fs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fs <- est_out(fs, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fs[, control_for_starting_distance := "yes"]

fx <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
     (1 | Year) + (1 | weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (scale(Covid)| IDLocality) + (1|sp_loc),
data = d[Country == "Finland"], REML = FALSE
)
est_fx <- est_out(fx, "Finland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_fx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hs <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hs <- est_out(hs, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hs[, control_for_starting_distance := "yes"]

hx <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Hungary"], REML = FALSE
)
est_hx <- est_out(hx, "Hungary: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_hx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
as <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
     (1 | Year) +(1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], REML = FALSE
)
est_as <- est_out(as, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_as[, control_for_starting_distance := "yes"]

ax <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_ax <- est_out(ax, "Australia: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(scale(Covid)|IDLocality)+(1|sp_loc)")
est_ax[, control_for_starting_distance := "no"]

# PL 
ps <- lmer(scale(log(FID)) ~
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_ps <- est_out(ps, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_ps[, control_for_starting_distance := "yes"]

px <- lmer(scale(log(FID)) ~
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(Covid) +
    (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = d[Country == "Poland"], REML = FALSE
)
est_px <- est_out(px, "Poland: (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)")
est_px[, control_for_starting_distance := "no"]

  # combine
    est_ms[, Country := 'All\n(mixed model)']
    est_mx[, Country := "All\n(mixed model)"]
    est_as[, Country := "Australia"]
    est_ax[, Country := "Australia"]
    est_cs[, Country := "Czechia"]
    est_cx[, Country := "Czechia"]
    est_hs[, Country := "Hungary"]
    est_hx[, Country := "Hungary"]
    est_ps[, Country := "Poland"]
    est_px[, Country := "Poland"]
    est_fs[, Country := "Finland"]
    est_fx[, Country := "Finland"]

    o = rbind(est_ms, est_mx, 
            est_as, est_ax, 
            est_cs, est_cx, 
            est_hs, est_hx,
            est_ps, est_px, 
            est_fs, est_fx)
    save(o, file = here::here('Data/dat_est_rev.Rdata'))
 
  #AIC(hs, hx)
  #AIC(cs, cx)
  #0AIC(ps, px)
  #AIC(as, ax)
  #AIC(fs, fx)

# supplementary Sable S_fid_period  
  ms_out <- m_out(name = "Table S1 - full a", dep = "Escape distance", model = ms, nsim = 5000)
  mx_out <- m_out(name = "Table S1 - full b", dep = "Escape distance", model = mx, nsim = 5000)
  cs_out <- m_out(name = "Table S1 - CZ a", dep = "Escape distancey", model = cs, nsim = 5000)
  cx_out <- m_out(name = "Table S1 - CZ b", dep = "Escape distance", model = cx, nsim = 5000)
  fs_out <- m_out(name = "Table S1 - FI a", dep = "Escape distance", model = fs, nsim = 5000)
  fx_out <- m_out(name = "Table S1 - FI b", dep = "Escape distance", model = fx, nsim = 5000)
  hs_out <- m_out(name = "Table S1 - HU a", dep = "Escape distance", model = hs, nsim = 5000)
  hx_out <- m_out(name = "Table S1 - HU b", dep = "Escape distance", model = hx, nsim = 5000)
  as_out <- m_out(name = "Table S1 - AU a", dep = "Escape distance", model = as, nsim = 5000)
  ax_out <- m_out(name = "Table S1 - AU b", dep = "Escape distance", model = ax, nsim = 5000)
  ps_out <- m_out(name = "Table S1 - PL a", dep = "Escape distance", model = ps, nsim = 5000)
  px_out <- m_out(name = "Table S1 - PL b", dep = "Escape distance", model = px, nsim = 5000)
  
  out_FID_c <- rbind(fs_out, fx_out, ps_out, px_out, cs_out, cx_out, hs_out, hx_out, as_out, ax_out, fill = TRUE)
  out_FID_c[is.na(out_FID_c)] <- ""
  out_FID_c$R2_mar = out_FID_c$R2_con = NULL
  out_FID_c[, effect := gsub("scale\\(Covid\\)", "Period", effect)]
  out_FID_c[, effect := gsub("scale\\(Year\\)", "year", effect)]
  out_FID_c[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
  out_FID_c[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
  out_FID_c[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
  out_FID_c[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
  out_FID_c[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
  out_FID_c[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
  out_FID_c[, effect := gsub("Species", "species", effect)]
  out_FID_c[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_c[, effect := gsub("IDLocality", "site", effect)]
  out_FID_c[, effect := gsub("sp_loc", "species within site", effect)]
  out_FID_c[, effect := gsub("site Period", "Period (slope) | site", effect)]
  fwrite(file = here::here("Outputs/Table_S1.csv"), out_FID_c)
load(here::here("Data/dat_est_rev.Rdata"))
o[predictor %in% c("scale(Covid)"), predictor := "Period"]
oo <- o[predictor %in% c("Period")]
oo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
  oo_s = oo[control_for_starting_distance == 'yes']
  met = summary(meta.summaries(d = oo_s$estimate, se = oo_s$sd, method = "fixed", weights = oo_s$N))$summci
  oo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

  oo_sx = oo[control_for_starting_distance == "no"]
  metx = summary(meta.summaries(d = oo_sx$estimate, se = oo_sx$sd, method = "fixed", weights = oo_sx$N))$summci
  oo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
  
  oo = rbind(oo, oo_met, oo_metx)
    
oo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

# prepare for adding N
oo[control_for_starting_distance == "no" | is.na(N), N := ""]
oo[, n_pos := 1.1]

width_ <- .5 # spacing between error bars

#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)

#p = 
ggplot(oo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
    geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
    geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
    # geom_point(position = ggstance::position_dodgev(.6)) +
    geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
    # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
    # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
    # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
    # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
    geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
    scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
    #scale_color_jama(guide = "none")+ #, palette = 'light'
    scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
    scale_x_continuous(breaks = round(seq(-0.6, 1.2, by = 0.3), 1)) +
    ylab("") +
    xlab("Standardized effect size of Period\n[on flight initiation distance]") +
    # coord_cartesian(xlim = c(-.15, .15)) +
    # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
    theme_bw() +
    theme(
        legend.position = "right",
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        # legend.spacing.y = unit(0.1, 'cm'),
        legend.key.height = unit(0.5, "line"),
        legend.margin = margin(0, 0, 0, 0),
        # legend.position=c(0.5,1.6),
        plot.title = element_text(color = "grey", size = 7),
        plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = ax_lines, size = 0.25),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
        # axis.text.x = element_text()
        axis.ticks.length = unit(1, "pt"),
        axis.text.x = element_text(, size = 6),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 7)
    )

ggsave(here::here("Outputs/Fig_rev_width_CustomLocusZoom_v2.png"), width = 8, height = 6.5, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_CustomJAMAv1.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_Okabe_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_UChicago_v3.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_d3_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_nejm_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jama_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_jco_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)
#ggsave("Outputs/Fig_rev_width_npg_v2.png", width = 8, height = 6, unit = "cm", dpi = 600)

# show used colors
#gg <- ggplot_build(p)
#col_ = unique(gg$data[[3]]["colour"])$colour
#show_col(col_)

Figure 1 | Change in avian tolerance towards humans before vs during the COVID-19 shutdowns. The dots with horizontal lines represent estimated standardised effect size and their 95% confidence intervals, the numbers sample sizes. For the country-specific and "All“, the effect sizes and 95% confidence intervals come from the joint posterior distribution of 5000 simulated values generated by the sim function from the arm package (Gelman et al., 2022) using the mixed model outputs controlled for starting distance of the observer (filled circles) or not (empty circles; Tables S1 and S2a). The models were further controlled for flock size, body size, temperature (also a proxy for a day within the breeding season: r Pearson = 0.48; Fig. S2), and time of a day, as well as for the non-independence of data points by fitting random intercepts of year, weekday, genus, species, species at a given day and year, country (in All - a global mixed model), site, and species within a site, while fitting Period as random slope within site (i.e. allowing for different Period effect at each site) and in All also within country. Fitting Period as random slope at other random intercepts produces similar results (see Fig. S1a). The multicollinearity was small as correlations between predictors were weak (Fig. S2). For the “Combined“, the estimate and 95% confidence interval represent the meta-analytical mean based on the country estimates and their standard deviation (from the country-specific models), and sample size per country. Note that effect sizes are small and estimates centre around zero.

Table S1 | Escape distance in relations to Period, given country

out_FID_c$error_structure = out_FID_c$response = NULL
out_FID_c[model!="", model:=c('Finland', 'Finland, without starting distance', 
                              'Poland', 'Poland, without starting distance',
                              'Czechia', 'Czechia, without starting distance',
                              'Hungary', 'Hungary, without starting distance',
                              'Australia', 'Australia, without starting distance')]
setnames(out_FID_c, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_FID_c %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model N type effect estimate lower upper
Finland 1019 fixed (Intercept) -0.035 -0.321 0.231
fixed starting distance (ln) 0.282 0.228 0.336
fixed flock size (ln) -0.049 -0.101 0.001
fixed body mass (ln) 0.091 -0.066 0.24
fixed time (sine of radians) 0.085 0.022 0.147
fixed time (cosine of radians) 0.031 -0.028 0.088
fixed temperaturre -0.079 -0.148 -0.009
fixed Periiod -0.15 -0.337 0.031
random species within day & year (Intercept) 10% 9% 11%
random species within site (Intercept) 6% 5% 7%
random site (Intercept) 1% 0% 3%
random site Periiod 1% 0% 3%
random species (Intercept) 4% 3% 5%
random genus (Intercept) 10% 7% 14%
random weekday (Intercept) 0% 0% 1%
random Year (Intercept) 5% 1% 11%
random Residual 61% 50% 72%
Finland, without starting distance 1019 fixed (Intercept) -0.125 -0.468 0.221
fixed flock size (ln) -0.018 -0.07 0.034
fixed body mass (ln) 0.2 0.03 0.375
fixed time (sine of radians) 0.093 0.026 0.161
fixed time (cosine of radians) 0.029 -0.032 0.089
fixed temperaturre -0.076 -0.148 -0.005
fixed Periiod -0.132 -0.366 0.101
random species within day & year (Intercept) 9% 7% 10%
random species within site (Intercept) 7% 6% 7%
random site (Intercept) 1% 1% 2%
random site Periiod 1% 1% 2%
random species (Intercept) 7% 5% 8%
random genus (Intercept) 10% 7% 13%
random weekday (Intercept) 1% 0% 1%
random Year (Intercept) 9% 2% 16%
random Residual 56% 44% 68%
Poland 762 fixed (Intercept) 0.145 -0.115 0.404
fixed starting distance (ln) 0.479 0.419 0.538
fixed flock size (ln) 0.012 -0.042 0.066
fixed body mass (ln) -0.07 -0.172 0.031
fixed time (sine of radians) -0.035 -0.117 0.05
fixed time (cosine of radians) -0.012 -0.101 0.078
fixed temperaturre -0.111 -0.182 -0.04
fixed Periiod 0.119 -0.112 0.353
random species within day & year (Intercept) 8% 8% 8%
random species (Intercept) 15% 12% 17%
random genus (Intercept) 0% 0% 0%
random weekday (Intercept) 0% 0% 1%
random Year (Intercept) 11% 7% 17%
random Residual 66% 57% 72%
Poland, without starting distance 762 fixed (Intercept) 0.246 -0.075 0.558
fixed flock size (ln) 0.054 -0.009 0.117
fixed body mass (ln) 0.046 -0.097 0.19
fixed time (sine of radians) -0.014 -0.112 0.077
fixed time (cosine of radians) -0.001 -0.104 0.099
fixed temperaturre -0.104 -0.187 -0.02
fixed Periiod 0.245 0.001 0.498
random species within day & year (Intercept) 11% 10% 11%
random species (Intercept) 26% 22% 29%
random genus (Intercept) 1% 1% 2%
random weekday (Intercept) 0% 0% 0%
random Year (Intercept) 8% 5% 14%
random Residual 53% 46% 60%
Czechia 2013 fixed (Intercept) 0.111 -0.257 0.481
fixed starting distance (ln) 0.392 0.351 0.437
fixed flock size (ln) 0.005 -0.03 0.037
fixed body mass (ln) 0.17 0.038 0.308
fixed time (sine of radians) 0.05 0.002 0.099
fixed time (cosine of radians) 0.054 0.009 0.097
fixed temperaturre 0.043 -0.008 0.093
fixed Periiod 0.036 -0.294 0.361
random species within day & year (Intercept) 3% 2% 3%
random species within site (Intercept) 3% 2% 3%
random species (Intercept) 13% 12% 13%
random site (Intercept) 2% 0% 3%
random site Periiod 2% 0% 3%
random genus (Intercept) 12% 11% 11%
random weekday (Intercept) 0% 0% 0%
random Year (Intercept) 15% 1% 33%
random Residual 52% 34% 70%
Czechia, without starting distance 2013 fixed (Intercept) -0.043 -0.969 0.867
fixed flock size (ln) 0.005 -0.031 0.04
fixed body mass (ln) 0.231 0.078 0.382
fixed time (sine of radians) 0.14 0.086 0.195
fixed time (cosine of radians) 0.123 0.075 0.172
fixed temperaturre 0.035 -0.027 0.098
fixed Periiod 0.223 -0.642 1.1
random species within day & year (Intercept) 2% 1% 5%
random species within site (Intercept) 1% 1% 2%
random species (Intercept) 6% 4% 10%
random site (Intercept) 1% -1% 1%
random site Periiod 1% -1% 1%
random genus (Intercept) 8% 5% 11%
random weekday (Intercept) 1% 0% 1%
random Year (Intercept) 54% 17% 73%
random Residual 26% 13% 57%
Hungary 1456 fixed (Intercept) 0.097 -0.174 0.384
fixed starting distance (ln) 0.485 0.438 0.531
fixed flock size (ln) -0.012 -0.051 0.031
fixed body mass (ln) -0.065 -0.184 0.059
fixed time (sine of radians) -0.077 -0.133 -0.019
fixed time (cosine of radians) -0.011 -0.06 0.037
fixed temperaturre -0.023 -0.083 0.037
fixed Periiod -0.091 -0.271 0.097
random species within day & year (Intercept) 5% 4% 5%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 3% 2% 3%
random genus (Intercept) 7% 5% 8%
random site (Intercept) 2% 0% 8%
random site Periiod 2% 0% 8%
random weekday (Intercept) 1% 0% 1%
random Year (Intercept) 5% 2% 8%
random Residual 76% 59% 84%
Hungary, without starting distance 1456 fixed (Intercept) 0.078 -0.336 0.494
fixed flock size (ln) -0.011 -0.057 0.035
fixed body mass (ln) 0.187 0.006 0.376
fixed time (sine of radians) -0.111 -0.18 -0.044
fixed time (cosine of radians) -0.007 -0.06 0.05
fixed temperaturre -0.025 -0.095 0.043
fixed Periiod -0.105 -0.307 0.109
random species within day & year (Intercept) 7% 5% 8%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 2% 2% 2%
random genus (Intercept) 18% 13% 18%
random site (Intercept) 1% 0% 12%
random site Periiod 1% 0% 12%
random weekday (Intercept) 1% 0% 1%
random Year (Intercept) 5% 3% 7%
random Residual 63% 42% 74%
Australia 1119 fixed (Intercept) -0.064 -0.202 0.081
fixed starting distance (ln) 0.502 0.452 0.553
fixed flock size (ln) 0.02 -0.026 0.068
fixed body mass (ln) -0.066 -0.163 0.029
fixed time (sine of radians) -0.032 -0.099 0.036
fixed time (cosine of radians) -0.026 -0.093 0.041
fixed temperaturre 0.009 -0.045 0.063
fixed Periiod -0.009 -0.08 0.062
random species within day & year (Intercept) 7% 6% 7%
random species within site (Intercept) 13% 11% 13%
random species (Intercept) 14% 11% 14%
random genus (Intercept) 2% 1% 2%
random site (Intercept) 0% -1% 8%
random site Periiod 0% -1% 8%
random weekday (Intercept) 0% 0% 0%
random Year (Intercept) 0% 0% 0%
random Residual 64% 51% 68%
Australia, without starting distance 1119 fixed (Intercept) -0.18 -0.432 0.058
fixed flock size (ln) 0.051 -0.002 0.104
fixed body mass (ln) 0.055 -0.073 0.184
fixed time (sine of radians) -0.017 -0.098 0.063
fixed time (cosine of radians) -0.006 -0.08 0.069
fixed temperaturre 0.004 -0.058 0.067
fixed Periiod 0.156 -0.007 0.331
random species within day & year (Intercept) 8% 5% 10%
random species within site (Intercept) 12% 8% 15%
random species (Intercept) 17% 14% 17%
random genus (Intercept) 5% 4% 5%
random site (Intercept) 1% -10% 15%
random site Periiod 1% -10% 15%
random weekday (Intercept) 0% 0% 0%
random Year (Intercept) 1% 0% 2%
random Residual 55% 36% 73%

Continuous variables were z-transformed. For model descriptions see legend of Fig. 1, for descriptions of variables Methods of the paper. Residual indicates residual variance.


Alternative models give similar results

  # prepare estimates Period
    # 01a all data, main text
      m1a <- lmer(scale(log(FID)) ~
        scale(log(SD)) +
        scale(log(FlockSize)) +
        scale(log(BodyMass)) +
        scale(sin(rad)) + scale(cos(rad)) +
        # scale(Day)+
        scale(Temp) +
        scale(Covid) +
        (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d, REML = FALSE,
       control = lmerControl(
         optimizer = "optimx", optCtrl = list(method = "nlminb")
       )
      )
     est_m1a = est_out(m1a, "01a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc)")
    # 01b all data, all random slopes - singularity 
      m1b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1 | weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE
                )
                  # # (Covid|IDLocality) +
      est_m1b = est_out(m1b, "01b) (1|Year) + (1|weekday) + (scale(Covid)|genus) + (scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01c all data, all random slopes, but some without cor to avoid singularity
      m1c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  #(1|Year) +(0+Covid|genus)+(0+Covid|Species)+(1|sp_day_year) + (Covid|Country) + (0+Covid|IDLocality) +(Covid|sp_loc)
                  (1|Year)+ (1|weekday) + (0+scale(Covid)|genus)+(0+scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE) # (0+scale(Covid)|genus) is zero 
      est_m1c = est_out(m1c, "01c) (1|Year) + (1|weekday) + (0+scale(Covid)|genus) + (0+scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01d all data, random slopes that allow for non-singular fit 
      m1d=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1|weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  data = d, REML = FALSE)# (Covid|IDLocality) +
      #d[,res := resid(mf)]
      est_m1d = est_out(m1d, "01d) (1|Year) + (1|weekday)+ (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc)")
    # 01e all data, random slopes that allow for non-singular fit with simple structure
      m1e=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  data = d, REML = FALSE)# (Covid|IDLocality) +
      est_m1e = est_out(m1e, "01e) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc)")

    # 02a) before & during > 4/species - main text 
      dx = dd[N_during>4 & N_before >4]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m2a <- lmer(scale(log(FID)) ~
        scale(log(SD)) +
        scale(log(FlockSize)) +
        scale(log(BodyMass)) +
        scale(sin(rad)) + scale(cos(rad)) +
        # scale(Day)+
        scale(Temp) +
        scale(Covid) +
        (1 | Year) + (1| weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d[Species %in% dx$Species], REML = FALSE,
       control = lmerControl(
         optimizer = "optimx", optCtrl = list(method = "nlminb")
       )
      )
      est_m2a = est_out(m2a, "02a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >4/species/period")
    
    # 02b) before & during > 4/species - singularity
      dx = dd[N_during > 4 & N_before > 4]
      m2b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m2b = est_out(m2b, "02b) (1|Year) + (1|weekday)+ (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >4/species/period")
    # 02c) before & during > 4/species, random slopes that allow for non-singular fit and simple structure
      dx = dd[N_during>4 & N_before >4]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m2c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year)+ (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m2c = est_out(m2c, "02c) (1|Year) + (1|weekday)+ (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >4/species/period")
    # 03a) before & during > 9/species - main text
      dx = dd[N_during > 9 & N_before > 9]
      m3a <- lmer(scale(log(FID)) ~
       scale(log(SD)) +
       scale(log(FlockSize)) +
       scale(log(BodyMass)) +
       scale(sin(rad)) + scale(cos(rad)) +
       # scale(Day)+
       scale(Temp) +
       scale(Covid) +
       (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(Covid) | Country) + (scale(Covid) | IDLocality) + (1 | sp_loc),
      data = d[Species %in% dx$Species], REML = FALSE,
      )
     est_m3a = est_out(m3a, "03a) (1|Year) + (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (1|sp_loc); >9/species/period")
    # 03b) before & during > 9/species - singularity 
      dx = dd[N_during>9 & N_before >9]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m3b=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday) + (scale(Covid)|genus)+(scale(Covid)|Species)+(1|sp_day_year) + (scale(Covid)|Country) +(scale(Covid)|IDLocality) + (scale(Covid)|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m3b = est_out(m3b, "03b) (1|Year) + (1|weekday) + (scale(Covid)|genus)+(scale(Covid)|Species) + (1|sp_day_year) + (scale(Covid)|Country) + (scale(Covid)|IDLocality) + (scale(Covid)|sp_loc); >9/species/period")
    # 03c) before & during > 9/species, random slopes that allow for non-singular fit and simple structure
      dx = dd[N_during>9 & N_before >9]
      #dxx = d[Species %in% dx$Species]
      #dx2 = dd[N_during>9 & N_before >9]
      m3c=lmer(scale(log(FID))~
                  scale(log(SD))+
                  scale(log(FlockSize))+
                  scale(log(BodyMass))+
                  scale(sin(rad)) + scale(cos(rad)) + 
                  #scale(Day)+
                  scale(Temp)+
                  scale(Covid)+
                  (1|Year) + (1| weekday) +(1|genus)+(1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc),
                  #(1|Year) + (1|genus)+(1|Species)+(1|sp_day_year) + (Covid|Country) + (Covid|sp_loc),
                  data = d[Species %in% dx$Species],
                  REML = FALSE) # (Covid|IDLocality) +
      est_m3c = est_out(m3c, '03c) (1|Year) + (1|genus) + (1|Species)+(1|sp_day_year) + (scale(Covid)|Country) + (1|IDLocality) + (1|sp_loc); >9/species/period')  
  
  # prepare estimates Stringency 
    m01a <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(StringencyIndex) +
      (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
      data = s, REML = FALSE)
    est_m01a = est_out(m01a, "01a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")

    m01b=lmer(scale(log(FID))~ 
        scale(Year)+ 
        scale(log(SD))+ 
        scale(log(FlockSize))+ 
        scale(log(BodyMass))+ 
        scale(sin(rad)) + scale(cos(rad)) +  
        #scale(Day)+ 
        scale(Temp)+ 
        scale(StringencyIndex)+ 
        (1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
        data = s, REML = FALSE
        )  
        # (1|Year) explains nothing - could stay 
     est_m01b = est_out(m01b, "01b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
     
     m01c=lmer(scale(log(FID))~
          scale(Year)+       
          scale(log(SD))+
          scale(log(FlockSize))+
          scale(log(BodyMass))+
          scale(sin(rad)) + scale(cos(rad)) + 
          #scale(Day)+
          scale(Temp)+
          scale(StringencyIndex)+
          (scale(StringencyIndex)|Country) + (1|IDLocality),  
          data = s, REML = FALSE) 
          # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
     est_m01c = est_out(m01c, "01c) (scale(StringencyIndex)|Country) + (1|IDLocality)")

     m02a <- lmer(scale(log(FID)) ~
       scale(Year) +
       scale(log(SD)) +
       scale(log(FlockSize)) +
       scale(log(BodyMass)) +
       scale(sin(rad)) + scale(cos(rad)) +
       # scale(Day)+
       scale(Temp) +
       scale(StringencyIndex) +
       (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
     data = s[Nsp>4], REML = FALSE
     )

     est_m02a = est_out(m02a, "02a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
    m02b=lmer(scale(log(FID))~ 
        scale(Year)+ 
        scale(log(SD))+ 
        scale(log(FlockSize))+ 
        scale(log(BodyMass))+ 
        scale(sin(rad)) + scale(cos(rad)) +  
        #scale(Day)+ 
        scale(Temp)+ 
        scale(StringencyIndex)+ 
        (1|weekday) + (scale(StringencyIndex)|genus)+(1|Species)+(1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc),
        data = s[Nsp>4], REML = FALSE,  
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
        )  
        # (1|Year) explains nothing - could stay 
     est_m02b = est_out(m02b, "02b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >4/species")
     
     m02c=lmer(scale(log(FID))~
          scale(Year)+       
          scale(log(SD))+
          scale(log(FlockSize))+
          scale(log(BodyMass))+
          scale(sin(rad)) + scale(cos(rad)) + 
          #scale(Day)+
          scale(Temp)+
          scale(StringencyIndex)+
          (scale(StringencyIndex)|Country) + (1|IDLocality),  
          data = s[Nsp>4], REML = FALSE) 
          # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
      est_m02c = est_out(m02c, "02c) (scale(StringencyIndex)|Country) + (1|IDLocality); >4/species")

      m03a <- lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
         data = s[Nsp > 9], REML = FALSE
         )

         est_m03a = est_out(m03a, "03a) (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")
         m03b = lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (1 | weekday) + (scale(StringencyIndex) | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
         data = s[Nsp > 9], REML = FALSE
         )
         est_m03b = est_out(m03b, "03b) (1|weekday) + (scale(StringencyIndex)|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc); >9/species")

        m03c = lmer(scale(log(FID)) ~
           scale(Year) +
           scale(log(SD)) +
           scale(log(FlockSize)) +
           scale(log(BodyMass)) +
           scale(sin(rad)) + scale(cos(rad)) +
           # scale(Day)+
           scale(Temp) +
           scale(StringencyIndex) +
           (scale(StringencyIndex) | Country) + (1 | IDLocality),
         data = s[Nsp > 9], REML = FALSE
         )
         # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
         est_m03c = est_out(m03c, "03c) (scale(StringencyIndex)|Country) + 1|IDLocality); >9/species")

  # prepare estimates google mobility 
    mg01a=lmer(scale(log(FID))~
      scale(Year)+       
      scale(log(SD))+
      scale(log(FlockSize))+
      scale(log(BodyMass))+
      scale(sin(rad)) + scale(cos(rad)) + 
      #scale(Day)+
      scale(Temp)+
      scale(parks_percent_change_from_baseline)+
      (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc),  
      data = ss, REML = FALSE
      ) 
      # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg01a = est_out(mg01a, "01a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc)")
    
    mg01b <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(parks_percent_change_from_baseline) +
     (1|weekday) + (scale(parks_percent_change_from_baseline)| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
     data = ss, REML = FALSE
     ) 

     est_mg01b = est_out(mg01b, "01b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc)")
     
    mg01c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss, REML = FALSE
    )
    est_mg01c = est_out(mg01c, "01c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)")

    mg02a = lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp>4], REML = FALSE
    )
    # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg02a = est_out(mg02a, "02a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >4/specie")

    mg02b <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp>4], REML = FALSE
    )

    est_mg02b = est_out(mg02b, "02b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >4/specie")

    mg02c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss[Nsp>4], REML = FALSE
    )
    est_mg02c = est_out(mg02c, "02c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >4/specie")

    mg03a = lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp > 9], REML = FALSE
    )
    # (1|Year), (1|genus), (1|sp_day_year), (1|sp_loc),  explain nothing - could stay
    est_mg03a = est_out(mg03a, "03a) (1|weekday) + (1|genus) +(1|Species) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality)+(1|sp_loc); >9/specie")

    mg03b <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (1 | weekday) + (scale(parks_percent_change_from_baseline) | genus) + (1 | Species) + (1 | sp_day_year) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
    data = ss[Nsp > 9], REML = FALSE
    )

    est_mg03b = est_out(mg03b, "03b) (1|weeekday) + (scale(parks_percent_change_from_baseline)|genus)+(1|Species)+(1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc); >9/specie")

    mg03c <- lmer(scale(log(FID)) ~
      scale(Year) +
      scale(log(SD)) +
      scale(log(FlockSize)) +
      scale(log(BodyMass)) +
      scale(sin(rad)) + scale(cos(rad)) +
      # scale(Day)+
      scale(Temp) +
      scale(parks_percent_change_from_baseline) +
      (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality),
    data = ss[Nsp > 9], REML = FALSE
    )
    est_mg03c = est_out(mg03c, "03c) (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality); >9/specie")

    # export
      save(file = here::here("Data/Fig_S2_estimates.Rdata"), 
      est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c, 
      est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c, 
      est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
     
     # prepare plot for Period
     xs = rbind(est_m1a, est_m1b, est_m1c, est_m1d, est_m1e, est_m2a, est_m2b, est_m2c, est_m3a, est_m3b, est_m3c)
     #gsub("scale\\(Covid\\)", "Period", "(scale(Covid)|Country)")
     xs[, model := gsub("scale\\(Covid\\)", "Period", xs$model)]
     xs[, model := gsub("Year", "year", xs$model)]
     xs[, model := gsub("Species", "species", xs$model)]
     xs[, model := gsub("sp_day_year", "species within day & year", xs$model)]
     xs[, model := gsub("IDLocality", "site", xs$model)]
     xs[, model := gsub("sp_loc", "species within site", xs$model)]

     gs2_ =
       ggplot(xs[predictor == "scale(Covid)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Period (before vs during shutdowns)", tag = "a)") + # title = "a)",Effect of Period (before/during shutdown)
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         #plot.title.position = "plot",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #gs2
     # ggsave(here::here('Outputs/Figure_Sy.png'),g, width = 30, height =5, units = 'cm')
     # prepare plot for Stringency
     xs0 = rbind(est_m01a, est_m01b, est_m01c, est_m02a, est_m02b, est_m02c, est_m03a, est_m03b, est_m03c)
     xs0[, model := gsub("scale\\(StringencyIndex\\)", "Stringency Index", model)]
     xs0[, model := gsub("Year", "year", model)]
     xs0[, model := gsub("Species", "species", model)]
     xs0[, model := gsub("sp_day_year", "species within day & year", model)]
     xs0[, model := gsub("IDLocality", "site", model)]
     xs0[, model := gsub("sp_loc", "species within site", model)]
     g0_ =
       ggplot(xs0[predictor == "scale(StringencyIndex)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Stringency index", tag = "b)") + # title = "b) Effect of ") +
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #g0
     # ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
     # prepare plot for Google
     xg0 = rbind(est_mg01a, est_mg01b, est_mg01c, est_mg02a, est_mg02b, est_mg02c, est_mg03a, est_mg03b, est_mg03c)
     xg0[, model := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", model)]
     xg0[, model := gsub("Year", "year", model)]
     xg0[, model := gsub("Species", "species", model)]
     xg0[, model := gsub("sp_day_year", "species within day & year", model)]
     xg0[, model := gsub("IDLocality", "site", model)]
     xg0[, model := gsub("sp_loc", "species within site", model)]
     gg0_ =
       ggplot(xg0[predictor == "scale(parks_percent_change_from_baseline)"], aes(y = model, x = estimate, col = model)) +
       geom_vline(xintercept = 0, col = "grey30", lty = 3) +
       geom_errorbar(aes(xmin = lwr, xmax = upr, col = model), width = 0, position = position_dodge(width = 0.01)) +
       # ggtitle ("Sim based")+
       geom_point(position = position_dodge(width = 0.01)) +
       # scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       # scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
       scale_color_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
       scale_y_discrete(limits = rev) +
       coord_fixed(ratio = 0.05, xlim = c(-0.23, 0.15)) +
       # scale_shape(guide = guide_legend(reverse = TRUE)) +
       # scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
       labs(y = NULL, x = "Google Mobility\n[Standardized effect sizes on\nflight initiation distances]", tag = "c)") + # title = "b) Effect of ") +
       # ylim(c(0,100))+
       # coord_flip()+
       theme_bw() +
       theme(
         legend.position = "none",
         plot.subtitle = element_text(size = 7),
         plot.title = element_text(size = 7),
         plot.tag = element_text(size = 7),
         legend.title = element_text(size = 7),
         legend.text = element_text(size = 6),
         ## legend.spacing.y = unit(0.1, 'cm'),
         legend.key.height = unit(0.5, "line"),
         # plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
         panel.grid = element_blank(),
         panel.border = element_blank(),
         panel.background = element_blank(),
         axis.line = element_line(colour = ax_lines, size = 0.25),
         axis.line.y = element_blank(),
         axis.ticks.y = element_blank(),
         axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
         axis.ticks.length = unit(1, "pt"),
         axis.text.x = element_text(colour = "black", size = 6),
         axis.text.y = element_text(colour = "black", size = 7),
         axis.title = element_text(size = 7)
       )
     #gg0
     # ggsave(here::here('Outputs/Figure_Sz.png'),g0, width = 30, height =5, units = 'cm')
     # combine
     grid.draw(rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)))

     if(save_plot==TRUE){
     ggsave(here::here("Outputs/Fig_S2_rev_v6.png"), rbind(ggplotGrob(gs2_), ggplotGrob(g0_), ggplotGrob(gg0_)), width = 30, height = 15, units = "cm")
     }

Figure S1 | Comparing estimates from alternative models. Changes in avian tolerance towards humans in response to (a) Period (before vs during the COVID-19 shutdowns) (b) stringency of governmental measures and (c) Google Mobility. The dots with horizontal lines represent the estimated standardised effect size and their 95% confidence intervals based on the joint posterior distribution of 5,000 simulated values generated by the sim function in R (Gelman et al. 2016) from the output of the mixed models (for details see Table S2). The name of each effect size highlights the corresponding model in Table S2a for (a), Table S2b for (b) and Table S2c for (c), the random structure of the specific model, if applicable, the condition used to reduce the dataset, and sample size. Depicted are effect sizes based on full (01) and reduced datasets with ≥5 (02) or ≥10 observations per species and period (03). For plots of model assumptions see https://doi.org/10.17605/OSF.IO/WUZH7 (Bulla et al. 2022). Note that effect sizes are small and estimates centre around zero.

Table S2a | Alternative models on escape distance given Period

    m1a_ = m_out(name = "Table S2a - 1a", dep = "Escape distance", model = m1a, nsim = 5000)
    m1b_ = m_out(name = "Table S2a - 1b", dep = "Escape distance", model = m1b, nsim = 5000)
    m1c_ = m_out(name = "Table S2a - 1c", dep = "Escape distance", model = m1c, nsim = 5000)
    m1d_ = m_out(name = "Table S2a - 1d", dep = "Escape distance", model = m1d, nsim = 5000)
    m1e_ = m_out(name = "Table S2a - 1e", dep = "Escape distance", model = m1e, nsim = 5000)

    m2a_ = m_out(name = "Table S2a - 2a", dep = "Escape distance", model = m2a, nsim = 5000)
    m2b_ = m_out(name = "Table S2a - 2b", dep = "Escape distance", model = m2b, nsim = 5000)
    m2c_ = m_out(name = "Table S2a - 2c", dep = "Escape distance", model = m2c, nsim = 5000)
    m3a_ = m_out(name = "Table S2a - 3a", dep = "Escape distance", model = m3a, nsim = 5000)
    m3b_ = m_out(name = "Table S2a - 3b", dep = "Escape distance", model = m3b, nsim = 5000)
    m3c_ = m_out(name = "Table S2a - 3c", dep = "Escape distance", model = m3c, nsim = 5000)
     
    out1 = rbind(m1a_, m1b_, m1c_, m1d_, m1e_, m2a_, m2b_, m2c_, m3a_, m3b_, m3c_, fill = TRUE)
    out1[is.na(out1)] = ""
    out1[, effect := gsub("scale\\(Covid\\)", "Period", effect)]
    out1[, effect := gsub("scale\\(Year\\)", "year", effect)]
    out1[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
    out1[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
    out1[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
    out1[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
    out1[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
    out1[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
    out1[, effect := gsub("Species", "species", effect)]
    out1[, effect := gsub("Country", "country", effect)]
    out1[, effect := gsub("Year", "year", effect)]
    out1[, effect := gsub("sp_day_year", "species within day & year", effect)]
    out1[, effect := gsub("IDLocality", "site", effect)]
    out1[, effect := gsub("sp_loc", "species within site", effect)]
    out1[type == "random" & grepl("Period", effect, fixed = TRUE), effect := paste("Period (slope) |", gsub(" Period", "", effect))]

    out1$R2_mar=out1$R2_con=NULL
    fwrite(file = here::here("Outputs/Table_S2a.csv"), out1)

    out1$response = out1$error_structure = NULL
    out1[model != "", model := paste0('0',substring(model, 12))]
    setnames(out1, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
    out1 %>%
       kbl() %>%
       kable_paper("hover", full_width = F)
model N type effect estimate lower upper
0 1a 6369 fixed (Intercept) 0.082 -0.129 0.286
fixed starting distance (ln) 0.477 0.453 0.501
fixed flock size (ln) -0.008 -0.027 0.011
fixed body mass (ln) 0.048 -0.02 0.111
fixed time (sine of radians) 0.022 -0.006 0.05
fixed time (cosine of radians) 0.028 0.006 0.051
fixed temperaturre 0 -0.029 0.028
fixed Period -0.033 -0.198 0.132
random species within day & year (Intercept) 8% 6% 9%
random species within site (Intercept) 5% 4% 5%
random site (Intercept) 1% 0% 5%
random Period (slope) | site 1% 0% 5%
random species (Intercept) 9% 8% 8%
random genus (Intercept) 6% 6% 6%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -3% 5%
random Period (slope) | country 1% -3% 5%
random year (Intercept) 3% 1% 4%
random Residual 65% 51% 76%
0 1b 6369 fixed (Intercept) 0.068 -0.155 0.294
fixed starting distance (ln) 0.478 0.455 0.502
fixed flock size (ln) -0.007 -0.026 0.012
fixed body mass (ln) 0.048 -0.017 0.113
fixed time (sine of radians) 0.021 -0.006 0.049
fixed time (cosine of radians) 0.027 0.004 0.05
fixed temperaturre -0.002 -0.03 0.028
fixed Period -0.048 -0.204 0.113
random species within day & year (Intercept) 12% 7% 14%
random species within site (Intercept) 0% -1% 0%
random Period (slope) | species within site 0% -1% 0%
random site (Intercept) 2% 0% 5%
random Period (slope) | site 2% 0% 5%
random species (Intercept) 0% -1% 8%
random Period (slope) | species 0% -1% 8%
random genus (Intercept) 0% 0% 5%
random Period (slope) | genus 0% 0% 5%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -4% 6%
random Period (slope) | country 1% -4% 6%
random year (Intercept) 3% 2% 4%
random Residual 79% 42% 95%
0 1c 6369 fixed (Intercept) 0.015 -0.148 0.181
fixed starting distance (ln) 0.494 0.469 0.517
fixed flock size (ln) -0.018 -0.037 0
fixed body mass (ln) 0.024 -0.007 0.054
fixed time (sine of radians) 0.015 -0.014 0.043
fixed time (cosine of radians) 0.027 0.004 0.051
fixed temperaturre -0.004 -0.034 0.026
fixed Period -0.041 -0.212 0.136
random species within day & year (Intercept) 12% 7% 13%
random species within site (Intercept) 1% 0% 12%
random Period (slope) | species within site 1% 0% 12%
random site (Intercept) 1% 1% 6%
random Period (slope) | site 1% 1% 6%
random Period (slope) | species 1% 1% 1%
random Period (slope) | genus 0% 0% 0%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -2% 4%
random Period (slope) | country 1% -2% 4%
random year (Intercept) 3% 2% 4%
random Residual 77% 46% 87%
0 1d 6369 fixed (Intercept) 0.084 -0.128 0.292
fixed starting distance (ln) 0.476 0.452 0.5
fixed flock size (ln) -0.008 -0.027 0.011
fixed body mass (ln) 0.048 -0.017 0.116
fixed time (sine of radians) 0.022 -0.006 0.05
fixed time (cosine of radians) 0.028 0.004 0.051
fixed temperaturre 0 -0.029 0.029
fixed Period -0.035 -0.2 0.132
random species within day & year (Intercept) 8% 6% 9%
random species within site (Intercept) 0% 0% 4%
random Period (slope) | species within site 0% 0% 4%
random site (Intercept) 1% 0% 5%
random Period (slope) | site 1% 0% 5%
random species (Intercept) 10% 9% 10%
random genus (Intercept) 7% 6% 6%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -3% 5%
random Period (slope) | country 1% -3% 5%
random year (Intercept) 3% 2% 4%
random Residual 66% 48% 78%
0 1e 6369 fixed (Intercept) 0.075 -0.125 0.282
fixed starting distance (ln) 0.478 0.454 0.5
fixed flock size (ln) -0.008 -0.027 0.011
fixed body mass (ln) 0.048 -0.019 0.116
fixed time (sine of radians) 0.022 -0.007 0.049
fixed time (cosine of radians) 0.024 0.001 0.047
fixed temperaturre -0.006 -0.034 0.023
fixed Period -0.033 -0.203 0.135
random species within day & year (Intercept) 7% 6% 8%
random species within site (Intercept) 5% 4% 5%
random site (Intercept) 6% 6% 6%
random species (Intercept) 8% 8% 9%
random genus (Intercept) 6% 5% 6%
random weekday (Intercept) 0% 0% 1%
random country (Intercept) 2% -3% 5%
random Period (slope) | country 2% -3% 5%
random year (Intercept) 2% 1% 4%
random Residual 62% 52% 71%
0 2a 5261 fixed (Intercept) 0.159 -0.183 0.512
fixed starting distance (ln) 0.462 0.437 0.489
fixed flock size (ln) -0.01 -0.031 0.011
fixed body mass (ln) -0.014 -0.124 0.096
fixed time (sine of radians) 0.028 -0.001 0.057
fixed time (cosine of radians) 0.036 0.012 0.06
fixed temperaturre -0.006 -0.038 0.027
fixed Period -0.036 -0.211 0.135
random species within day & year (Intercept) 7% 4% 9%
random species within site (Intercept) 3% 2% 4%
random site (Intercept) 1% -1% 2%
random Period (slope) | site 1% -1% 2%
random species (Intercept) 10% 8% 9%
random genus (Intercept) 5% 4% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -6% 20%
random Period (slope) | country 2% -6% 20%
random year (Intercept) 3% 2% 3%
random Residual 64% 34% 84%
0 2b 5261 fixed (Intercept) 0.16 -0.192 0.493
fixed starting distance (ln) 0.462 0.436 0.489
fixed flock size (ln) -0.01 -0.031 0.01
fixed body mass (ln) -0.011 -0.12 0.1
fixed time (sine of radians) 0.027 -0.002 0.056
fixed time (cosine of radians) 0.036 0.011 0.061
fixed temperaturre -0.005 -0.037 0.026
fixed Period -0.035 -0.211 0.137
random species within day & year (Intercept) 9% 4% 11%
random species within site (Intercept) 0% 0% 2%
random Period (slope) | species within site 0% 0% 2%
random site (Intercept) 1% -1% 2%
random Period (slope) | site 1% -1% 2%
random species (Intercept) 0% -1% 5%
random Period (slope) | species 0% -1% 5%
random genus (Intercept) 0% 0% 6%
random Period (slope) | genus 0% 0% 6%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 3% -6% 17%
random Period (slope) | country 3% -6% 17%
random year (Intercept) 4% 2% 3%
random Residual 78% 30% 100%
0 2c 5261 fixed (Intercept) 0.152 -0.183 0.503
fixed starting distance (ln) 0.464 0.438 0.491
fixed flock size (ln) -0.01 -0.03 0.01
fixed body mass (ln) -0.008 -0.121 0.099
fixed time (sine of radians) 0.026 -0.002 0.056
fixed time (cosine of radians) 0.032 0.008 0.057
fixed temperaturre -0.012 -0.043 0.02
fixed Period -0.032 -0.217 0.145
random species within day & year (Intercept) 8% 5% 9%
random species within site (Intercept) 3% 2% 4%
random site (Intercept) 4% 3% 4%
random species (Intercept) 8% 7% 7%
random genus (Intercept) 6% 4% 6%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -6% 17%
random Period (slope) | country 2% -6% 17%
random year (Intercept) 3% 2% 4%
random Residual 63% 39% 81%
0 3a 5107 fixed (Intercept) 0.133 -0.216 0.476
fixed starting distance (ln) 0.462 0.435 0.489
fixed flock size (ln) -0.01 -0.03 0.011
fixed body mass (ln) 0.006 -0.104 0.113
fixed time (sine of radians) 0.027 -0.003 0.057
fixed time (cosine of radians) 0.034 0.009 0.06
fixed temperaturre -0.007 -0.04 0.026
fixed Period -0.028 -0.196 0.144
random species within day & year (Intercept) 7% 4% 9%
random species within site (Intercept) 3% 2% 4%
random site (Intercept) 1% -1% 2%
random Period (slope) | site 1% -1% 2%
random species (Intercept) 21% 16% 18%
random genus (Intercept) 0% 0% 0%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -5% 19%
random Period (slope) | country 2% -5% 19%
random year (Intercept) 3% 2% 3%
random Residual 60% 32% 79%
0 3b 5107 fixed (Intercept) 0.141 -0.197 0.491
fixed starting distance (ln) 0.461 0.433 0.489
fixed flock size (ln) -0.01 -0.03 0.011
fixed body mass (ln) 0.013 -0.104 0.13
fixed time (sine of radians) 0.027 -0.003 0.058
fixed time (cosine of radians) 0.034 0.009 0.059
fixed temperaturre -0.006 -0.04 0.026
fixed Period -0.03 -0.203 0.141
random species within day & year (Intercept) 9% 3% 11%
random species within site (Intercept) 0% 0% 2%
random Period (slope) | species within site 0% 0% 2%
random site (Intercept) 1% -1% 2%
random Period (slope) | site 1% -1% 2%
random species (Intercept) 0% -1% 9%
random Period (slope) | species 0% -1% 9%
random genus (Intercept) 0% 0% 3%
random Period (slope) | genus 0% 0% 3%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 3% -5% 16%
random Period (slope) | country 3% -5% 16%
random year (Intercept) 4% 2% 3%
random Residual 79% 28% 100%
0 3c 5107 fixed (Intercept) 0.127 -0.221 0.478
fixed starting distance (ln) 0.464 0.437 0.49
fixed flock size (ln) -0.01 -0.031 0.01
fixed body mass (ln) 0.005 -0.105 0.12
fixed time (sine of radians) 0.027 -0.003 0.056
fixed time (cosine of radians) 0.028 0.004 0.053
fixed temperaturre -0.014 -0.046 0.018
fixed Period -0.031 -0.208 0.145
random species within day & year (Intercept) 7% 5% 9%
random species within site (Intercept) 3% 2% 4%
random site (Intercept) 5% 3% 5%
random species (Intercept) 15% 11% 15%
random genus (Intercept) 0% 0% 0%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -6% 18%
random Period (slope) | country 2% -6% 18%
random year (Intercept) 3% 2% 4%
random Residual 62% 39% 78%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance. Note that (1a) model is the one reported in the main text.

Table S2b | Alternative models on escape distance given Stringency

    m01a_ = m_out(name = "Table S2b - 1a", dep = "Escape distance", model = m01a, nsim = 5000)
    m01b_ = m_out(name = "Table S2b - 1b", dep = "Escape distance", model = m01b, nsim = 5000)
    m01c_ = m_out(name = "Table S2b - 1c", dep = "Escape distance", model = m01c, nsim = 5000)

    m02a_ = m_out(name = "Table S2b - 2a", dep = "Escape distance", model = m02a, nsim = 5000)
    m02b_ = m_out(name = "Table S2b - 2b", dep = "Escape distance", model = m02b, nsim = 5000)
    m02c_ = m_out(name = "Table S2b - 2c", dep = "Escape distance", model = m02c, nsim = 5000)

    m03a_ = m_out(name = "Table S2b - 3a", dep = "Escape distance", model = m03a, nsim = 5000)
    m03b_ = m_out(name = "Table S2b - 3b", dep = "Escape distance", model = m03b, nsim = 5000)
    m03c_ = m_out(name = "Table S2b - 3c", dep = "Escape distance", model = m03c, nsim = 5000)

    out2 = rbind(m01a_, m01b_, m01c_, m02a_, m02b_, m02c_, m03a_, m03b_, m03c_, fill = TRUE)
    out2[is.na(out2)] = ""
    out2[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
    out2[, effect := gsub("scale\\(Year\\)", "year", effect)]
    out2[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
    out2[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
    out2[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
    out2[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
    out2[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
    out2[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
    out2[, effect := gsub("Species", "species", effect)]
    out2[, effect := gsub("Country", "country", effect)]
    out2[, effect := gsub("Year", "year", effect)]
    out2[, effect := gsub("sp_day_year", "species within day & year", effect)]
    out2[, effect := gsub("IDLocality", "site", effect)]
    out2[, effect := gsub("sp_loc", "species within site", effect)]
    out2[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
    out2$R2_mar = out2$R2_con = NULL

    fwrite(file = here::here("Outputs/Table_S2b.csv"), out2)

    out2$response = out2$error_structure = NULL
    out2[model != "", model := paste0("0", substring(model, 12))]
    setnames(out2, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
    out2 %>%
      kbl() %>%
      kable_paper("hover", full_width = F)
model N type effect estimate lower upper
0 1a 3676 fixed (Intercept) -0.074 -0.345 0.173
fixed year 0.02 -0.028 0.068
fixed starting distance (ln) 0.492 0.46 0.524
fixed flock size (ln) -0.003 -0.03 0.023
fixed body mass (ln) 0.021 -0.046 0.091
fixed time (sine of radians) -0.02 -0.061 0.022
fixed time (cosine of radians) -0.005 -0.04 0.032
fixed temperaturre -0.007 -0.049 0.036
fixed stringency index 0.018 -0.157 0.182
random species within day & year (Intercept) 5% 3% 7%
random species within site (Intercept) 5% 3% 6%
random species (Intercept) 7% 6% 9%
random site (Intercept) 8% 7% 9%
random genus (Intercept) 4% 4% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -16% 16%
random stringency index (slope) | country 2% -16% 16%
random Residual 67% 45% 96%
0 1b 3676 fixed (Intercept) -0.073 -0.319 0.173
fixed year 0.019 -0.028 0.065
fixed starting distance (ln) 0.493 0.462 0.525
fixed flock size (ln) -0.003 -0.029 0.023
fixed body mass (ln) 0.017 -0.052 0.085
fixed time (sine of radians) -0.02 -0.061 0.021
fixed time (cosine of radians) -0.004 -0.039 0.03
fixed temperaturre -0.008 -0.049 0.033
fixed stringency index 0.006 -0.16 0.174
random species within day & year (Intercept) 5% 3% 7%
random species within site (Intercept) 5% 3% 7%
random species (Intercept) 7% 5% 8%
random site (Intercept) 8% 6% 10%
random genus (Intercept) 0% -1% 4%
random stringency index (slope) | genus 0% -1% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -16% 15%
random stringency index (slope) | country 2% -16% 15%
random Residual 70% 44% 100%
0 1c 3676 fixed (Intercept) -0.103 -0.278 0.066
fixed year -0.003 -0.05 0.044
fixed starting distance (ln) 0.549 0.518 0.581
fixed flock size (ln) -0.023 -0.05 0.003
fixed body mass (ln) -0.024 -0.053 0.005
fixed time (sine of radians) -0.039 -0.081 0.002
fixed time (cosine of radians) 0.005 -0.031 0.043
fixed temperaturre 0.003 -0.039 0.045
fixed stringency index 0.021 -0.14 0.181
random site (Intercept) 10% 10% 10%
random country (Intercept) 1% -7% 8%
random stringency index (slope) | country 1% -7% 8%
random Residual 88% 74% 103%
0 2a 3573 fixed (Intercept) -0.154 -0.455 0.133
fixed year 0.022 -0.025 0.071
fixed starting distance (ln) 0.483 0.451 0.516
fixed flock size (ln) 0.003 -0.023 0.029
fixed body mass (ln) -0.024 -0.1 0.052
fixed time (sine of radians) -0.02 -0.061 0.022
fixed time (cosine of radians) -0.008 -0.043 0.029
fixed temperaturre -0.008 -0.05 0.034
fixed stringency index 0.024 -0.153 0.191
random species within day & year (Intercept) 5% 3% 7%
random species within site (Intercept) 4% 3% 7%
random site (Intercept) 7% 5% 9%
random species (Intercept) 11% 8% 13%
random genus (Intercept) 2% 2% 2%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 3% -20% 20%
random stringency index (slope) | country 3% -20% 20%
random Residual 65% 39% 103%
0 2b 3573 fixed (Intercept) -0.151 -0.441 0.146
fixed year 0.021 -0.029 0.07
fixed starting distance (ln) 0.484 0.452 0.517
fixed flock size (ln) 0.003 -0.023 0.03
fixed body mass (ln) -0.03 -0.107 0.044
fixed time (sine of radians) -0.02 -0.06 0.021
fixed time (cosine of radians) -0.007 -0.043 0.03
fixed temperaturre -0.008 -0.051 0.034
fixed stringency index 0.015 -0.159 0.189
random species within day & year (Intercept) 5% 3% 8%
random species within site (Intercept) 4% 3% 7%
random site (Intercept) 7% 5% 9%
random species (Intercept) 11% 8% 13%
random genus (Intercept) 0% -1% 2%
random stringency index (slope) | genus 0% -1% 2%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 3% -21% 19%
random stringency index (slope) | country 3% -21% 19%
random Residual 67% 39% 106%
0 2c 3573 fixed (Intercept) -0.109 -0.286 0.064
fixed year -0.001 -0.049 0.048
fixed starting distance (ln) 0.543 0.511 0.576
fixed flock size (ln) -0.018 -0.045 0.009
fixed body mass (ln) -0.027 -0.055 0.003
fixed time (sine of radians) -0.038 -0.08 0.003
fixed time (cosine of radians) 0.002 -0.036 0.039
fixed temperaturre 0.002 -0.04 0.047
fixed stringency index 0.023 -0.136 0.19
random site (Intercept) 10% 10% 10%
random country (Intercept) 1% -8% 8%
random stringency index (slope) | country 1% -8% 8%
random Residual 88% 74% 105%
0 3a 3425 fixed (Intercept) -0.137 -0.459 0.173
fixed year 0.026 -0.023 0.077
fixed starting distance (ln) 0.48 0.448 0.512
fixed flock size (ln) 0.001 -0.027 0.028
fixed body mass (ln) 0.023 -0.061 0.11
fixed time (sine of radians) -0.012 -0.054 0.03
fixed time (cosine of radians) -0.009 -0.046 0.027
fixed temperaturre -0.011 -0.054 0.032
fixed stringency index 0.02 -0.153 0.187
random species within day & year (Intercept) 5% 3% 7%
random species within site (Intercept) 5% 3% 7%
random site (Intercept) 8% 5% 10%
random species (Intercept) 9% 7% 10%
random genus (Intercept) 3% 3% 3%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 2% -21% 20%
random stringency index (slope) | country 2% -21% 20%
random Residual 67% 39% 105%
0 3b 3425 fixed (Intercept) -0.14 -0.456 0.181
fixed year 0.027 -0.022 0.075
fixed starting distance (ln) 0.482 0.449 0.514
fixed flock size (ln) 0.001 -0.026 0.029
fixed body mass (ln) 0.019 -0.067 0.106
fixed time (sine of radians) -0.011 -0.053 0.032
fixed time (cosine of radians) -0.009 -0.045 0.028
fixed temperaturre -0.012 -0.055 0.033
fixed stringency index 0.013 -0.155 0.189
random species within day & year (Intercept) 5% 3% 7%
random species within site (Intercept) 4% 3% 7%
random site (Intercept) 8% 5% 11%
random species (Intercept) 10% 7% 11%
random genus (Intercept) 0% 0% 2%
random stringency index (slope) | genus 0% 0% 2%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 3% -23% 21%
random stringency index (slope) | country 3% -23% 21%
random Residual 67% 36% 110%
0 3c 3425 fixed (Intercept) -0.107 -0.288 0.076
fixed year 0.005 -0.045 0.054
fixed starting distance (ln) 0.541 0.508 0.574
fixed flock size (ln) -0.02 -0.047 0.007
fixed body mass (ln) -0.02 -0.048 0.01
fixed time (sine of radians) -0.03 -0.073 0.012
fixed time (cosine of radians) -0.001 -0.038 0.037
fixed temperaturre -0.002 -0.045 0.042
fixed stringency index 0.023 -0.137 0.183
random site (Intercept) 11% 11% 11%
random country (Intercept) 1% -8% 8%
random stringency index (slope) | country 1% -8% 8%
random Residual 87% 73% 105%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance. Note that (1a) model is the one reported in the main text.

Table S2c | Alternative models on escape distance given Google Mobility

    # google
     mg01a_ = m_out(name = "Table S2c - 1a", dep = "Escape distance", model = mg01a, nsim = 5000)
     mg01b_ = m_out(name = "Table S2c - 1b", dep = "Escape distance", model = mg01b, nsim = 5000)
     mg01c_ = m_out(name = "Table S2c - 1c", dep = "Escape distance", model = mg01c, nsim = 5000)

     mg02a_ = m_out(name = "Table S2c - 2a", dep = "Escape distance", model = mg02a, nsim = 5000)
     mg02b_ = m_out(name = "Table S2c - 2b", dep = "Escape distance", model = mg02b, nsim = 5000)
     mg02c_ = m_out(name = "Table S2c - 2c", dep = "Escape distance", model = mg02c, nsim = 5000)

     mg03a_ = m_out(name = "Table S2c - 3a", dep = "Escape distance", model = mg03a, nsim = 5000)
     mg03b_ = m_out(name = "Table S2c - 3b", dep = "Escape distance", model = mg03b, nsim = 5000)
     mg03c_ = m_out(name = "Table S2c - 3c", dep = "Escape distancey", model = mg03c, nsim = 5000)

     out3 = rbind(mg01a_, mg01b_, mg01c_, mg02a_, mg02b_, mg02c_, mg03a_, mg03b_, mg03c_, fill = TRUE)
     out3[is.na(out3)] = ""
      out3[, effect := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", effect)]
      out3[, effect := gsub("scale\\(Year\\)", "year", effect)]
      out3[, effect := gsub("scale\\(log\\(SD\\)\\)", "starting distance (ln)", effect)]
      out3[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
      out3[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
      out3[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
      out3[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
      out3[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
       out3[, effect := gsub("Species", "species", effect)]
       out3[, effect := gsub("Country", "country", effect)]
       out3[, effect := gsub("Year", "year", effect)]
       out3[, effect := gsub("sp_day_year", "species within day & year", effect)]
       out3[, effect := gsub("IDLocality", "site", effect)]
       out3[, effect := gsub("sp_loc", "species within site", effect)]
       out3[type == "random" & grepl("Google Mobility", effect, fixed = TRUE), effect := paste("Google Mobility (slope) |", gsub(" Google Mobility", "", effect))]
       out3$R2_mar = out3$R2_con = NULL
     fwrite(file = here::here("Outputs/Table_S2c.csv"), out3)
     out3$response = out3$error_structure = NULL
     out3[model != "", model := paste0("0", substring(model, 12))]
     setnames(out3, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
     out3 %>%
          kbl() %>%
          kable_paper("hover", full_width = F)
model N type effect estimate lower upper
0 1a 3644 fixed (Intercept) 0.025 -0.097 0.146
fixed year 0.019 -0.023 0.06
fixed starting distance (ln) 0.499 0.468 0.531
fixed flock size (ln) -0.001 -0.028 0.026
fixed body mass (ln) 0.025 -0.042 0.092
fixed time (sine of radians) -0.009 -0.05 0.033
fixed time (cosine of radians) 0.01 -0.023 0.043
fixed temperaturre 0.021 -0.017 0.059
fixed Google Mobility -0.1 -0.227 0.026
random species within day & year (Intercept) 6% 5% 6%
random species within site (Intercept) 5% 4% 5%
random species (Intercept) 7% 7% 8%
random site (Intercept) 9% 8% 10%
random genus (Intercept) 3% 3% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -3% 4%
random Google Mobility (slope) | country 1% -3% 4%
random Residual 68% 61% 78%
0 1b 3644 fixed (Intercept) 0.024 -0.093 0.141
fixed year 0.019 -0.021 0.061
fixed starting distance (ln) 0.498 0.466 0.53
fixed flock size (ln) -0.001 -0.027 0.025
fixed body mass (ln) 0.029 -0.038 0.095
fixed time (sine of radians) -0.01 -0.05 0.032
fixed time (cosine of radians) 0.01 -0.023 0.042
fixed temperaturre 0.022 -0.017 0.06
fixed Google Mobility -0.101 -0.232 0.026
random species within day & year (Intercept) 6% 5% 6%
random species within site (Intercept) 5% 4% 5%
random species (Intercept) 7% 7% 7%
random site (Intercept) 9% 9% 9%
random genus (Intercept) 0% -1% 4%
random Google Mobility (slope) | genus 0% -1% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -3% 4%
random Google Mobility (slope) | country 1% -3% 4%
random Residual 71% 59% 81%
0 1c 3644 fixed (Intercept) -0.057 -0.143 0.026
fixed year 0.013 -0.027 0.053
fixed starting distance (ln) 0.559 0.528 0.59
fixed flock size (ln) -0.021 -0.048 0.004
fixed body mass (ln) -0.023 -0.052 0.005
fixed time (sine of radians) -0.033 -0.075 0.009
fixed time (cosine of radians) 0.02 -0.014 0.056
fixed temperaturre 0.031 -0.007 0.069
fixed Google Mobility -0.087 -0.206 0.034
random site (Intercept) 11% 10% 12%
random country (Intercept) 0% -1% 2%
random Google Mobility (slope) | country 0% -1% 2%
random Residual 88% 83% 92%
0 2a 3545 fixed (Intercept) -0.027 -0.171 0.11
fixed year 0.019 -0.025 0.061
fixed starting distance (ln) 0.491 0.458 0.522
fixed flock size (ln) 0.004 -0.023 0.031
fixed body mass (ln) -0.019 -0.094 0.054
fixed time (sine of radians) -0.01 -0.051 0.032
fixed time (cosine of radians) 0.009 -0.025 0.043
fixed temperaturre 0.019 -0.018 0.059
fixed Google Mobility -0.098 -0.232 0.034
random species within day & year (Intercept) 6% 5% 6%
random species within site (Intercept) 4% 4% 5%
random site (Intercept) 9% 8% 9%
random species (Intercept) 10% 9% 11%
random genus (Intercept) 2% 1% 2%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 0% -5% 5%
random Google Mobility (slope) | country 0% -5% 5%
random Residual 69% 59% 80%
0 2b 3545 fixed (Intercept) -0.035 -0.168 0.098
fixed year 0.019 -0.024 0.061
fixed starting distance (ln) 0.49 0.458 0.522
fixed flock size (ln) 0.004 -0.022 0.03
fixed body mass (ln) -0.03 -0.101 0.043
fixed time (sine of radians) -0.011 -0.052 0.031
fixed time (cosine of radians) 0.009 -0.025 0.044
fixed temperaturre 0.021 -0.018 0.059
fixed Google Mobility -0.09 -0.222 0.045
random species within day & year (Intercept) 5% 5% 6%
random species within site (Intercept) 4% 4% 5%
random site (Intercept) 8% 8% 8%
random species (Intercept) 16% 14% 17%
random genus (Intercept) 0% 0% 0%
random Google Mobility (slope) | genus 0% 0% 0%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 0% -5% 4%
random Google Mobility (slope) | country 0% -5% 4%
random Residual 66% 56% 77%
0 2c 3545 fixed (Intercept) -0.056 -0.141 0.028
fixed year 0.014 -0.024 0.053
fixed starting distance (ln) 0.551 0.519 0.583
fixed flock size (ln) -0.017 -0.044 0.011
fixed body mass (ln) -0.026 -0.055 0.004
fixed time (sine of radians) -0.033 -0.074 0.009
fixed time (cosine of radians) 0.018 -0.018 0.053
fixed temperaturre 0.031 -0.008 0.068
fixed Google Mobility -0.084 -0.209 0.032
random site (Intercept) 11% 9% 12%
random country (Intercept) 0% -1% 2%
random Google Mobility (slope) | country 0% -1% 2%
random Residual 89% 84% 93%
0 3a 3399 fixed (Intercept) -0.024 -0.178 0.131
fixed year 0.027 -0.017 0.071
fixed starting distance (ln) 0.487 0.454 0.52
fixed flock size (ln) 0.002 -0.025 0.029
fixed body mass (ln) 0.025 -0.061 0.113
fixed time (sine of radians) 0 -0.041 0.045
fixed time (cosine of radians) 0.007 -0.028 0.043
fixed temperaturre 0.022 -0.018 0.06
fixed Google Mobility -0.11 -0.235 0.019
random species within day & year (Intercept) 6% 5% 6%
random species within site (Intercept) 4% 4% 5%
random site (Intercept) 9% 8% 9%
random species (Intercept) 7% 6% 9%
random genus (Intercept) 3% 3% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 0% -5% 4%
random Google Mobility (slope) | country 0% -5% 4%
random Residual 70% 60% 81%
0 3b 3399 fixed (Intercept) -0.011 -0.135 0.112
fixed year 0.031 -0.012 0.073
fixed starting distance (ln) 0.487 0.454 0.518
fixed flock size (ln) 0.003 -0.024 0.03
fixed body mass (ln) 0.036 -0.048 0.122
fixed time (sine of radians) -0.002 -0.044 0.042
fixed time (cosine of radians) 0.007 -0.026 0.042
fixed temperaturre 0.023 -0.017 0.063
fixed Google Mobility -0.128 -0.258 -0.003
random species within day & year (Intercept) 6% 5% 6%
random species within site (Intercept) 4% 4% 5%
random site (Intercept) 9% 9% 9%
random species (Intercept) 8% 7% 9%
random genus (Intercept) 0% -1% 4%
random Google Mobility (slope) | genus 0% -1% 4%
random weekday (Intercept) 0% 0% 0%
random country (Intercept) 1% -2% 3%
random Google Mobility (slope) | country 1% -2% 3%
random Residual 69% 58% 79%
0 3c 3399 fixed (Intercept) -0.068 -0.155 0.02
fixed year 0.023 -0.019 0.063
fixed starting distance (ln) 0.549 0.517 0.583
fixed flock size (ln) -0.019 -0.047 0.009
fixed body mass (ln) -0.019 -0.049 0.011
fixed time (sine of radians) -0.024 -0.066 0.018
fixed time (cosine of radians) 0.015 -0.021 0.053
fixed temperaturre 0.03 -0.01 0.069
fixed Google Mobility -0.091 -0.204 0.023
random site (Intercept) 11% 10% 12%
random country (Intercept) 0% 0% 2%
random Google Mobility (slope) | country 0% 0% 2%
random Residual 88% 83% 91%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance. Note that (1a) model is the one reported in the main text.

# TODO: modelAss

Correlations among predictors

d[, sin_rad := sin(rad)]
d[, cos_rad := cos(rad)]

dp <- d[, c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day")]
setnames(dp, old = c("SD_ln", "flock_ln", "body_ln", "sin_rad", "cos_rad", "Temp", "Day"), new = c("Starting distance\nln(m)", "Flock size\nln(m)", "Body mass\nln(m)", "Sine\nof radians", "Cosine\nof radians", "Temperature\n°C", "Day"))

#if (save_plot == TRUE) {
#  png(here::here("Outputs/Fig_S1_rev.png"), width = 19, height = 19, units = "cm", bg = "transparent", res = 600)
#  chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
#  mtext("Single observations", side = 3, line = 3)
#  dev.off()
#}
chart.Correlation(dp, histogram = TRUE, pch = 19, alpha = 0.5)
mtext("Single observations", side = 3, line = 3)

Figure S2 | Pairwise correlations among fixed effects used in this study. On the diagonal: histograms and density lines (red) for each variable. Above diagonal: Pearson’s correlation with stars indicating significance. Below diagonal: the bivariate scatterplots, with each dot representing a single observation and a red line representing smoothed fit. Created with chart.Correlation function from R-package PerformanceAnalytics (Peterson & Carl 2020).


Species-site specific distributions

px = pp[N_during > 4 & N_before > 4]
dxx = d[paste(IDLocality, Species) %in% paste(px$IDLocality, px$Species)]
# table(dxx$IDLocality, dxx$Year)
#length(unique(px$IDLocality))
#length(unique(px$Species))
#sum(px$N_during) + sum(px$N_before)

dxx[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
dxx[, sp_C_loc2 := paste(gsub("[_]", " ", Species), Country, IDLocality, sep = "\n")]
dxx[, genus := sub("_.*", "", Species)]
dxx[Covid == 0, period := "before COVID-19"]
dxx[Covid == 1, period := "during COVID-19"]

col3_ <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ <- col3_[3:7]

  ggplot(dxx, aes(x = as.factor(Year), y = FID, col = Country, fill = period)) +
  geom_boxplot(outlier.size = 0.5) +
  facet_wrap(~sp_C_loc2) +
  scale_y_continuous("Flight initiation distance [m]", trans = "log10") +
  scale_x_discrete("Year", guide = guide_axis(angle = 45)) +
  # scale_color_continuous() +
  scale_colour_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("white", "lightgrey")) +
  theme_MB +
  theme(
    plot.title = element_text(size = 7),
    strip.background = element_blank(),
    strip.text.x = element_text(size = 5, color = "grey30", margin = margin(1, 1, 1, 1, "mm")),
    # panel.spacing = unit(1, "mm"),
    legend.position = "none", # c(1, 0.01),
    legend.justification = c(1, 0),
    legend.title = element_blank(),
    # legend.spacing.y = unit(-0.78, "cm")
    # legend.spacing.y = unit(0.02, "cm") use if LOESS smooth text as legend
    legend.spacing.y = unit(-0.9, "cm"),
    axis.text.x = element_text(colour = "grey30", size = 6),
    axis.text.y = element_text(colour = "grey30", size = 6)
  )

if(save_plot == TRUE){ggsave(here::here("Outputs/Fig_2_rev_v2.png"), width = 18, height = 16, units = "cm")}

Figure 2 | Temporal variation in avian tolerance toward humans across species. Each heading denotes the scientific name of the species, country and unique site ID within each country. Boxplots outline colour highlights country (as in Fig. 1), fill colour indicates Period (white: before the COVID-19 shutdowns; grey: during the COVID-19 shutdowns). Boxplots depict median (horizontal line inside the box), the 25th and 75th percentiles (box) ± 1.5 times the interquartile range or the minimum/maximum value, whichever is smaller (bars), and the outliers (dots). Included are only species–site combinations with ≥5 observations per Period. Y-axis is on the log-scale. Note the lack of consistent shutdowns effects within and between species as well as within and between the countries.


Between/within-genus variation

# prepare data for fig 3 & s3
dxx <- d[paste(IDLocality, Species) %in% paste(pp$IDLocality, pp$Species)]
# length(dxx[, unique(paste(sp_loc, Covid))])
m <- lm(log(FID) ~ log(SD), dxx)
dxx[, resid_FID := resid(m)]
a <- dxx[, .(mean(resid_FID), sd(resid_FID), mean(FID), .N), by = .(Country, IDLocality, genus, Species, sp_loc, Covid)]
setnames(a, old = c("V1", "V2", "V3"), new = c("resid_FID_avg", "SD", "FID_avg"))
a[is.na(SD), SD := 0]

aw <- reshape(a, idvar = c("Country", "IDLocality", "genus", "Species", "sp_loc"), timevar = "Covid", direction = "wide")
aw[, Species := gsub("[_]", " ", Species)]
aw <- merge(aw, t, all.x = TRUE)
# table(aw$Family)

x <- aw[, .N, by = Species]
# x[order(Species)]
aw[, genus2 := genus]
aw[Species %in% x[N %in% c(1, 2), Species], genus2 := "other"]

# aw[genus2 == "Phoenicurus", unique(Species)]

ph[genus2 == "Motacilla" | uid %in% c("67a9ecfd-58ba-44a4-9986-243b6e610419"), uid := "cf522e02-35cc-44f5-841c-0e642987c2e4"]
ph[genus2 == "Sylvia", uid := "67a9ecfd-58ba-44a4-9986-243b6e610419"]

ph[, size := 0.2]
ph[genus2 %in% c("Anas", "Columba", "Dendrocopos", "Sturnus"), size := c(0.25, 0.25, 0.15, 0.1)]
ph[, FID_avg.0 := 1.5]
ph[, FID_avg.1 := 20]
ph[genus2 %in% c("Anas", "Columba"), FID_avg.0 := c(1.7, 1.7)]

ph[, resid_FID_avg.0 := -1.7]
ph[, resid_FID_avg.1 := 0.7]
ph[genus2 %in% c("Anas", "Columba"), resid_FID_avg.0 := c(-1.6, -1.6)]

ph[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

aw[, genus2 := factor(genus2, levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))]

# Fig 3 and left panel of S3- plot from files
anas <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Anas.png"))))
columba <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Columba.png"))))
Dendrocopos <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Dendrocopos.png"))))
Larus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Larus_flip.png"))))
Picus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Picus.png"))))
Motacilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Motacilla.png"))))
Erithacus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Erithacus.png"))))
Phoenicurus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Phoenicurus.png"))))
Turdus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Turdus.png"))))
Sylvia <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sylvia.png"))))
Parus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Parus_flip.png"))))
Sitta <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sitta.png"))))
Pica <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Pica.png"))))
Garrulus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Garrulus.png"))))
Corvus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Corvus.png"))))
Sturnus <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Sturnus.png"))))
Passer <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Passer_flip.png"))))
Fringilla <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/Fringilla.png"))))
other <- rasterGrob(change_col("#CCCCCC", readPNG(here::here("Data/Pics/other_flip.png"))))

ann_text <- data.frame(
  FID_avg.0 = 8, FID_avg.1 = 10, lab = "Text",
  genus2 = factor("Anas", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)
ann_text2 <- data.frame(
  FID_avg.0 = 6, FID_avg.1 = 3, lab = "Text",
  genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other"))
)

aw2 <- data.frame(FID_avg.0 = c(11.25, 11.25), FID_avg.1 = c(3.5, 5.8), genus2 = factor("Larus", levels = c("Anas", "Larus", "Columba", "Dendrocopos", "Picus", "Motacilla", "Erithacus", "Phoenicurus", "Turdus", "Sylvia", "Parus", "Sitta", "Pica", "Garrulus", "Corvus", "Sturnus", "Passer", "Fringilla", "other")))

col3_ <- c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ <- col3_[3:7]

aw[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]
g_gen <-
  ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) +
  # geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
  # geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
  # geom_point(pch = 21, alpha = 0.7, aes(col = Country)) +
  annotation_custom2(anas, data = ph[genus2 == "Anas"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
  annotation_custom2(Larus, data = ph[genus2 == "Larus"], xmin = 0.05, xmax = 0.5, ymax = 2.6) +
  annotation_custom2(columba, data = ph[genus2 == "Columba"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
  annotation_custom2(Dendrocopos, data = ph[genus2 == "Dendrocopos"], xmin = 0.05, xmax = 0.25, ymax = 2.6) +
  annotation_custom2(Picus, data = ph[genus2 == "Picus"], xmin = 0.05, xmax = 0.4, ymax = 2.7) +
  annotation_custom2(Motacilla, data = ph[genus2 == "Motacilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Erithacus, data = ph[genus2 == "Erithacus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
  annotation_custom2(Phoenicurus, data = ph[genus2 == "Phoenicurus"], xmin = 0.05, xmax = 0.35, ymax = 2.7) +
  annotation_custom2(Turdus, data = ph[genus2 == "Turdus"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Sylvia, data = ph[genus2 == "Sylvia"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Parus, data = ph[genus2 == "Parus"], xmin = 0.05, xmax = 0.42, ymax = 2.7) +
  annotation_custom2(Sitta, data = ph[genus2 == "Sitta"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(Pica, data = ph[genus2 == "Pica"], xmin = 0.05, xmax = 0.5, ymax = 2.5) +
  annotation_custom2(Garrulus, data = ph[genus2 == "Garrulus"], xmin = 0.05, xmax = 0.6, ymax = 2.7) +
  annotation_custom2(Corvus, data = ph[genus2 == "Corvus"], xmin = 0.05, xmax = 0.4, ymax = 2.55) +
  annotation_custom2(Sturnus, data = ph[genus2 == "Sturnus"], xmin = 0.05, xmax = 0.24, ymax = 2.65) +
  annotation_custom2(Passer, data = ph[genus2 == "Passer"], xmin = 0.05, xmax = 0.36, ymax = 2.65) +
  annotation_custom2(Fringilla, data = ph[genus2 == "Fringilla"], xmin = 0.05, xmax = 0.5, ymax = 2.7) +
  annotation_custom2(other, data = ph[genus2 == "other"], xmin = 0.05, xmax = 0.45, ymax = 2.4) +
  geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
  # ggtitle ("Sim based")+
  geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
  geom_line(data = aw2, col = "grey80", lwd = 0.25) +
  geom_text(data = ann_text, label = "No difference", col = "grey80", angle = 45, size = 2) +
  geom_text(data = ann_text2, label = "Species mean / site", col = "grey60", size = 2, ) +
  facet_wrap(~genus2) +
  # geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) + # ,
  # scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = FALSE)) +
  scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  # labs(title = "Species means per sampling location")+
  theme_MB +
  theme(
    plot.title = element_text(size = 7),
    strip.background = element_blank(),
    # panel.spacing = unit(1, "mm"),
    legend.position = c(1, 0.025),
    legend.justification = c(1, -0.05)
  )
gg_gen <- ggplotGrob(g_gen) # gg$layout$name
ggx_gen <- gtable_filter_remove(gg_gen,
  name = paste0("axis-b-", c(2, 4), "-4"),
  trim = FALSE
)
if (save_plot == TRUE) {
  ggsave(here::here("Outputs/Fig_3_width-122mm_col_grey_rev.png"), ggx_gen, width = 4.8, height = 4.5, dpi = 600) # 12.2cm # with label inside
}
grid.draw(ggx_gen)

Figure 3 | Avian tolerance towards humans before and during the COVID-19 shutdowns according to genera. Dots represent means or single escape distance observations of species at specific sites (e.g. park or cemetery) with data for both periods (before and during the shutdowns) and not corrected for other factors such as starting distance of the observer (plot corrected for starting distance gives similar patterns: Fig. S3). Dot colour highlights the country. Dotted lines indicate no difference; dots above the lines indicate lower tolerance towards humans (i.e. longer escape distances), dots below the lines indicate higher tolerance during than before the COVID-19 shutdowns. Panels are ordered according to evolutionary history of birds with top left panels representing the oldest genera, and bottom right, the youngest. Panel titled ‘other’ contains genera with only one or two data points. The axes are on the log-scale. For a species-specific figure, see Fig. S4. Silhouette of Garrulus glandarius, Motacilla alba, Picus viridis, Phoenicurus ochruros, Sylvia borin were drawn by Martin Bulla, Erithacus rubecula drawn by Rebecca Groom, and Fringilla coelebs and Sturnus vulgaris by Maxime Dahirel and all are available at PhyloPic under Creative Commons Attribution 3.0 Unported licence. The remaining silhouettes are available at PhyloPic under the Public Domain Dedication 1.0 license.

# Fig S3 right panel
g_2 <-
  ggplot(aw, aes(x = resid_FID_avg.0, y = resid_FID_avg.1)) +
  geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
  geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
  facet_wrap(~genus2) +
  # geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) +
  # scale_fill_viridis(discrete = TRUE, guide = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous("Before COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
  scale_y_continuous("During COVID-19 shutdown - residual escape distance", expand = c(0, 0)) +
  # labs(title = "Species means per sampling location")+
  theme_MB +
  theme(
    plot.title = element_text(size = 7),
    strip.background = element_blank(),
    legend.position = "none",
    # legend.position = c(1, 0),
    legend.justification = c(1, 0)
  )

g_g2 <- ggplotGrob(g_2) # gg$layout$name
g_gx2 <- gtable_filter_remove(g_g2,
  name = paste0("axis-b-", c(2, 4), "-4"),
  trim = FALSE
)
# Fig S5 combine
grid.draw(cbind(ggx_gen, g_gx2, size = "last"))

if (save_plot == TRUE) {
  ggsave(here::here("Outputs/Fig_S3_rev_v3.png"), cbind(ggx_gen, g_gx2, size = "last"), width = 4.8 * 2, height = 4.5, dpi = 600)
}

Figure S3 | Comparison of genus-specific flight initiation distance (left) and residual flight initiation distance (right) before and during the COVID-19 shutdown. Left panel is a copy of a main text Fig. 3 (see there for details). The escape distance represents the raw data that can be confounded by the observers starting distance (for our data the r Pearson = 0.58. Right panel dots depict residual escape distances from a model with flight initiation distance (ln-transformed) as a response and starting distance (ln-transformed) as a predictor, i.e. the dots represent before and during shutdowns values that are controlled for starting distance. Note that such control for starting distance (right) has little influence on the depicted relationships and whereas some genera might have decreased their escape distance, other likely increased it and yet most kept the escape distance similar across the two periods. Indeed, the Pearson’s correlation coefficient for escape distance (ln-scale) and residual escape distance was 0.8 for single values and 0.74 for the species means per sampling location.

Difference in number of observations per species and site before and during shutdowns:

ggplot(aw, aes(x = N.0 - N.1)) +
  geom_histogram() + xlab('Before minus during shutdowns\n[# observations')

# nrow(aw[abs(N.0 - N.1) > 2])
# nrow(aw[!abs(N.0 - N.1) > 2])

Between/within-species variation

aw[, sp2 := gsub(" ", "\n", Species)]
ann_text <- data.frame(
  FID_avg.0 = 8, FID_avg.1 = 10, lab = "Text",
  Species = factor("Aegithalos caudatus", levels = levels(as.factor(aw$Species)))
)
ann_text$sp2 = gsub(" ", "\n", ann_text$Species)
g3 <-
  ggplot(aw, aes(x = FID_avg.0, y = FID_avg.1)) +
  # geom_errorbar(aes(ymin = FID_avg.1-SD.1, ymax = FID_avg.1+SD.1, col = Country), width = 0) +
  # geom_errorbar(aes(xmin = FID_avg.0-SD.0, xmax = FID_avg.0+SD.0, col = Country), width = 0) +
  # geom_point(pch = 21, alpha = 0.7, aes(col = Country)) +
  geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = "grey20", alpha = 0.8) + # col = "white") +
  # ggtitle ("Sim based")+
  geom_abline(intercept = 0, slope = 1, lty = 3, col = "grey80") +
  geom_text(data = ann_text, label = "No difference", col = "grey80", angle = 45, size = 2) +
  facet_wrap(~sp2) +
  # geom_phylopic(data = o, aes(image = uid),  color = "grey80", size = o$size) + # ,
  # scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = FALSE))  +
  scale_fill_manual(values = col3__, guide = guide_legend(reverse = TRUE)) +
  scale_x_continuous("Before COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  scale_y_continuous("During COVID-19 shutdown - flight initiation distance [m]", expand = c(0, 0), trans = "log10") +
  labs(title = "Species means per sampling location") +
  theme_MB +
  theme(
    plot.title = element_text(size = 7),
    strip.background = element_blank(),
    # panel.spacing = unit(1, "mm"),
    legend.position = c(0.96, 0.0),
    legend.justification = c(1, 0)
  )
gg3 <- ggplotGrob(g3) # gg2$layout$name
ggx3 <- gtable_filter_remove(gg3,
  name = c(paste0("axis-b-", c(2, 4), "-7"), "axis-b-6-6"),
  trim = FALSE
)
grid.draw(ggx3)

if (save_plot == TRUE) {
  ggsave(here::here("Outputs/Fig_S4_species_rev_v4.png"), ggx3, width = 13.5, height = 17.5, unit = "cm", dpi = 600) # 11.43cm
}

Figure S4 | Avian tolerance towards humans before and during the COVID-19 shutdowns according to species. Dots represent means or single escape distance observations of species at specific sites (e.g. specific park or cemetery) with data for both periods (i.e. before and during the shutdowns) and not corrected for other factors such as starting distance of the observer. Dot colour highlights the country. Dotted lines indicate no difference, dots above the lines indicate lower tolerance towards humans (i.e. longer escape distances), dots below the lines indicate lower tolerance before than during the COVID-19 shutdowns. Panels are ordered alphabetically. The axes are on the log-scale.


Exploration of Google Mobility

g_ <- fread(here::here("Data/google_mobility.txt")) # fwrite(d, here::here('Data/data.txt'), sep ='\t')
g_[, Year := as.integer(substring(date, nchar(date) - 3, nchar(date)))]
g_[nchar(date) == 9, date := paste0("0", date)]
g_[, date_ := as.Date(date, format = "%d.%m.%Y")]
g_[, Day := yday(date_)]
setnames(g_, old = "country_region", new = "Country")
g_[, weekday := weekdays(date_)]

g0 = ggplot(g_, aes(x = parks_percent_change_from_baseline, fill = factor(Year))) +
  geom_histogram(position = "dodge") +
  # scale_y_continuous(trans = 'log')+
  scale_fill_manual(values = c("orange", "skyblue", "black"), guide = 'none') +
  geom_vline(xintercept = 0, lty = 3, col = "red") +
  labs(subtitle = "Distribution") +
  xlab( 'Google Mobility\n[% change in human presence]') +
  ylab( 'Count') +
  facet_wrap(~Country, nrow = 5)

g1 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
  geom_line() +
  facet_wrap(~Country, nrow = 5) +
  labs(subtitle = "Raw data", xlab = 'Day\n ', y = 'Google Mobility\n[% change in human presence]') +
  # scale_y_continuous(trans = 'log')+
  coord_cartesian(ylim = c(-100, 300))+
  scale_color_manual(values = c("orange", "skyblue", "black"), guide = 'none')

g2 = ggplot(g_, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
     stat_smooth() +
     facet_wrap(~Country, nrow = 5) +
     labs(subtitle = "Loess", xlab = 'Day\n ') +
     # scale_y_continuous(trans = 'log')+
     coord_cartesian(ylim = c(-100, 300))+
     scale_color_manual(values = c("orange", "skyblue", "black"), name = 'Year')+
      theme(
         axis.text.y = element_blank(),
         axis.title.y = element_blank()
       )
  ggarrange(
    g0, g1, g2,
    ncol = 3, widths = c(0.9, 1, 1.1)
  )

 ggsave(here::here("Outputs/Fig_4_rev.png"), width = 8*2.54, height = 5*2.54, unit = "cm", dpi = 600)

Figure 4 | Changes in human presence in parks within and between years and countries. Left plots represent distribution (histograms) of human presence (Google Mobility), middle plots the raw data, right plots locally estimated scatterplot smoothing. Note that Google Mobility data are not freely available for years before the COVID-19 pandemic (i.e. before 2020) but 2022 was a year without shutdowns in the studied countries. For weekday-specific pattern see Fig. S5 below.

g[, weekday := factor(weekday, levels = (c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")))]
ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
    geom_line() +
    facet_grid(rows = vars(Country), cols = vars(weekday)) +
    # scale_y_continuous(trans = 'log')+
    scale_color_manual(values = c("orange", "skyblue", "black"))

ggsave(here::here("Outputs/Fig_S5_rev.png"), width = 8 * 2.54, height = 6 * 2.54, unit = "cm", dpi = 600)

Figure S5 | Changes in human presence (Google Mobility) in parks across weekdays and years. Depicted are raw data connected by lines. Note that Google Mobility data were not freely available for years before the COVID-19 pandemic (i.e. before 2020) but 2022 was a year without shutdowns in the studied countries.


Google Mobility ~ stringency

# Predictions 
l = list()
 sc = s[Country == "Czechia"]
 cz <- lmer(parks_percent_change_from_baseline ~
    StringencyIndex + 
   (scale(StringencyIndex)|weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sc, REML = FALSE
 )
 bsim <- sim(cz, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(StringencyIndex = seq(min(sc$StringencyIndex), max(sc$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country = 'Czechia'
 l[[1]] = newD
 
 s[, year_weekday :=paste(Year, weekday)]
 sf = s[Country == "Finland"]
 fi <- lmer(parks_percent_change_from_baseline ~
    Year+
    StringencyIndex + 
   (scale(StringencyIndex)|year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sf, REML = FALSE
 )
 bsim <- sim(fi, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sf$Year), StringencyIndex = seq(min(sf$StringencyIndex), max(sf$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country = 'Finland'
 newD$Year = NULL
 l[[2]] = newD

 sh <- s[Country == "Hungary"]
 hu <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sh, REML = FALSE
 )

 bsim <- sim(hu, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sh$Year), StringencyIndex = seq(min(sh$StringencyIndex), max(sh$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Hungary"
 newD$Year <- NULL
 l[[3]] <- newD
 
 sp <- s[Country == "Poland"]
 pl <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sp, REML = FALSE
 )
 bsim <- sim(pl, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sp$Year), StringencyIndex = seq(min(sp$StringencyIndex), max(sp$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Poland"
 newD$Year <- NULL
 l[[4]] <- newD

 sa <- s[Country == "Australia"]
 au <- lmer(parks_percent_change_from_baseline ~
     Year +
     StringencyIndex +
     (scale(StringencyIndex) | year_weekday),
 # (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
 data = sa, REML = FALSE
 )

 bsim <- sim(au, n.sim = nsim)
 v <- apply(bsim@fixef, 2, quantile, prob = c(0.5))
 newD <- data.frame(Year = mean(sa$Year), StringencyIndex = seq(min(sa$StringencyIndex), max(sa$StringencyIndex), length.out = 100)) # values to predict for
 X <- model.matrix(~ Year + StringencyIndex, data = newD) # exactly the model which was used has to be specified here
 newD$pred <- (X %*% v)
 predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
 for (j in 1:nsim) predmatrix[, j] <- (X %*% bsim@fixef[j, ])
 newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
 newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
 newD$pred <- apply(predmatrix, 1, quantile, prob = 0.5)
 newD$Country <- "Australia"
 newD$Year <- NULL
 l[[5]] <- newD

# Figure G_S
g_s = data.table(do.call(rbind,l))
g_s[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia")))]

col3_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1]
col3__ = col3_[3:7]
#p = 
ggplot(g_s, aes(x = StringencyIndex, y = pred, col = Country)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = Country, color = NULL), alpha = .15) +
  geom_jitter(aes(y = parks_percent_change_from_baseline, fill = Country), data = s, pch = 21, col = 'grey20', width = 0.7, height = 3, alpha = 0.5) +
  geom_line(lwd = 1) +
  labs(subtitle = "Mixed model per country predicitons", y = "Google Mobiligy\n[% change from baseline]", x = "Stringency Index") +
   # scale_color_locuszoom()+
   # scale_fill_locuszoom(guide = "none")
  scale_x_continuous(breaks = round(seq(25, 75, by = 25), 1)) +
  scale_y_continuous(breaks = round(seq(-100, 200, by = 50), 1)) +
  #scale_y_continuous(breaks = round(seq(-100, 175, by = 25), 1)) +
  scale_colour_manual(values = col3__, guide = guide_legend(reverse = TRUE, override.aes = list(size = 0)), 
            labels = paste("<span style='color:",
                                   col3__,
                                   "'>",
                                   levels(sp$Country),
                                   "</span>")
            ) +
  scale_fill_manual(values = col3__, guide = "none") +
  theme_bw() +
  theme(
    legend.text = element_markdown(size = 6),
    #legend.position = "right",
    legend.title = element_blank(),
    # legend.spacing.y = unit(0.1, 'cm'),
    legend.key.height = unit(0.5, "line"),
    legend.key.size = unit(0, "line"),
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin=margin(-10,1,-10,-10),
    # legend.position=c(0.5,1.6),
    plot.title = element_text(color = "grey", size = 7),
    plot.subtitle = element_text(color = "grey60", size = 6),
    plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
    panel.grid = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = ax_lines, size = 0.25),
    axis.ticks = element_line(colour = ax_lines, size = 0.25),
    # axis.text.x = element_text()
    axis.ticks.length = unit(1, "pt"),
    axis.text = element_text(, size = 6),
    axis.title = element_text(size = 7)
  )

if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_5_rev.png"), width = 7, height = 6, unit = "cm", dpi = 600)
}

Figure 5 | Association between human presence in parks (Google Mobility) and stringency of antipandemic governmental restrictions (stringency index). Lines with shaded areas represent predicted relationships from country-specific mixed effect models controlled for the year and non-independence of data points by including weekday within the year as random intercept and stringency index as a random slope (Table S3). Dots represent raw data, jittered to increase visibility. Colours indicate country. Note the generally negative but weak association between human presence and stringency index.

Table S3 | Google Mobility in relation to stringency index

ll = list()
s[, year_weekday := paste(Year, weekday)]

sf = s[Country == "Finland"]
fi <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sf, REML = FALSE
)
ll[[1]] = m_out(name = "Table S3 - FI", dep = "Google Mobility", model = fi, nsim = 5000)

sp <- s[Country == "Poland"]
pl <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sp, REML = FALSE
)
ll[[2]] = m_out(name = "Table S3 - PL", dep = "Google Mobility", model = pl, nsim = 5000)

sc = s[Country == "Czechia"]
cz <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(StringencyIndex) +
  (scale(StringencyIndex) | weekday),
data = sc, REML = FALSE
)
ll[[3]] = m_out(name = "Table S3 - CZ", dep = "Google Mobility", model = cz, nsim = 5000)

sh <- s[Country == "Hungary"]
hu <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
data = sh, REML = FALSE
)
ll[[4]] = m_out(name = "Table S3 - HU", dep = "Google Mobility", model = hu, nsim = 5000)

sa <- s[Country == "Australia"]
au <- lmer(scale(parks_percent_change_from_baseline) ~
  scale(Year) +
  scale(StringencyIndex) +
  (scale(StringencyIndex) | year_weekday),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = sa, REML = FALSE
)
ll[[5]] = m_out(name = "Table S3 - AU", dep = "Google Mobility", model = au, nsim = 5000)

out_g_s = data.table(do.call(rbind, ll))
out_g_s[is.na(out_g_s)] <- ""
out_g_s$R2_mar = out_g_s$R2_con = NULL
out_g_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
out_g_s[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
out_g_s[, effect := gsub("year_weekday", "weekday within year", effect)]
out_g_s[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
fwrite(file = here::here("Outputs/Table_S3_rev.csv"), out_g_s)

out_g_s$error_structure = out_g_s$response = NULL
out_g_s[model!="", model:=c('Finland', 
                              'Poland', 
                              'Czechia', 
                              'Hungary', 
                              'Australia')]
setnames(out_g_s, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_g_s %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model N type effect estimate lower upper
Finland 322 fixed (Intercept) -0.27 -1.02 0.529
fixed year 0.936 0.441 1.44
fixed stringency index -0.261 -1.023 0.492
random weekday within year (Intercept) 22% 47% 54%
random stringency index (slope) | weekday within year 22% 47% 54%
random Residual 55% -7% 5%
Poland 329 fixed (Intercept) 0.023 -0.627 0.626
fixed year 0.76 0.015 1.531
fixed stringency index 0.016 -0.359 0.4
random weekday within year (Intercept) 47% 44% 49%
random stringency index (slope) | weekday within year 47% 44% 49%
random Residual 6% 2% 13%
Czechia 1054 fixed (Intercept) 0.009 -0.521 0.548
fixed stringency index -0.229 -0.642 0.191
random weekday (Intercept) 23% -182% 40%
random stringency index (slope) | weekday 23% -182% 40%
random Residual 54% 21% 464%
Hungary 874 fixed (Intercept) 0.249 -0.255 0.766
fixed year 0.002 -0.442 0.465
fixed stringency index -1.269 -1.969 -0.574
random weekday within year (Intercept) 48% 21% 49%
random stringency index (slope) | weekday within year 48% 21% 49%
random Residual 4% 2% 59%
Australia 1065 fixed (Intercept) 0.519 -0.011 1.053
fixed year 0.443 0.241 0.645
fixed stringency index -0.537 -0.848 -0.246
random weekday within year (Intercept) -31% 44% 65%
random stringency index (slope) | weekday within year -31% 44% 65%
random Residual 162% -30% 13%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance.


Effect sizes for stringency or Google Mobility

# predictions for fig and table for stringency
 # full
 mss <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(StringencyIndex) +
     (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1 | IDLocality) + (1 | sp_loc),
 data = s, REML = FALSE,
 control = lmerControl(
     optimizer = "optimx", optCtrl = list(method = "nlminb")
 )
 )
 est_mss <- est_out(mss, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
 est_mss[, control_for_starting_distance := "yes"]

 msx <- lmer(scale(log(FID)) ~
    scale(Year) +
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (scale(StringencyIndex) | Country) + (1| IDLocality) + (1 | sp_loc),
 data = s, REML = FALSE,
 control = lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_msx <- est_out(msx, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(StringencyIndex)|Country) + (1|IDLocality) +(1|sp_loc)")
est_msx[, control_for_starting_distance := "no"]


# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  css <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = s[Country == "Czechia"], REML = FALSE
    )
  est_css <- est_out(css, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_css[, control_for_starting_distance := "yes"]
  
  csx <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = s[Country == "Czechia"], REML = FALSE
    )
  est_csx <- est_out(csx, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_csx[, control_for_starting_distance := "no"]

# FI
fss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fss <- est_out(fss, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fss[, control_for_starting_distance := "yes"]

fsx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = s[Country == "Finland"], REML = FALSE
)
est_fsx <- est_out(fsx, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fsx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) +  (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hss <- est_out(hss, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hss[, control_for_starting_distance := "yes"]

hsx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Hungary"], REML = FALSE
)
est_hsx <- est_out(hsx, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hsx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ass <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
     (1|weekday) +  (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], REML = FALSE
)
est_ass <- est_out(ass, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ass[, control_for_starting_distance := "yes"]

asx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_asx <- est_out(asx, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_asx[, control_for_starting_distance := "no"]

# PL 
pss <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_pss <- est_out(pss, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pss[, control_for_starting_distance := "yes"]

psx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = s[Country == "Poland"], REML = FALSE
)
est_psx <- est_out(psx, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_psx[, control_for_starting_distance := "no"]

  # combine
    est_mss[, Country := 'All\n(mixed model)']
    est_msx[, Country := "All\n(mixed model)"]
    est_ass[, Country := "Australia"]
    est_asx[, Country := "Australia"]
    est_css[, Country := "Czechia"]
    est_csx[, Country := "Czechia"]
    est_hss[, Country := "Hungary"]
    est_hsx[, Country := "Hungary"]
    est_pss[, Country := "Poland"]
    est_psx[, Country := "Poland"]
    est_fss[, Country := "Finland"]
    est_fsx[, Country := "Finland"]

    os = rbind(est_mss, est_msx, 
            est_ass, est_asx, 
            est_css, est_csx, 
            est_hss, est_hsx,
            est_pss, est_psx, 
            est_fss, est_fsx)
    save(os, file = here::here('Data/dat_est_Stringency_rev.Rdata'))

# estimates for table
  mss_out <- m_out(name = "Table S4 - full a", dep = "Escape distance", model = mss, nsim = 5000)
  msx_out <- m_out(name = "Table S4 - full b", dep = "Escape distance", model = msx, nsim = 5000)
  css_out <- m_out(name = "Table S4 - CZ a", dep = "Escape distance", model = css, nsim = 5000)
  csx_out <- m_out(name = "Table S4 - CZ b", dep = "Escape distance", model = csx, nsim = 5000)
  fss_out <- m_out(name = "Table S4 - FI a", dep = "Escape distance", model = fss, nsim = 5000)
  fsx_out <- m_out(name = "Table S4 - FI b", dep = "Escape distance", model = fsx, nsim = 5000)
  hss_out <- m_out(name = "Table S4 - HU a", dep = "Escape distance", model = hss, nsim = 5000)
  hsx_out <- m_out(name = "Table S4 - HU b", dep = "Escape distance", model = hsx, nsim = 5000)
  ass_out <- m_out(name = "Table S4 - AU a", dep = "Escape distance", model = ass, nsim = 5000)
  asx_out <- m_out(name = "Table S4 - AU b", dep = "Escape distancey", model = asx, nsim = 5000)
  pss_out <- m_out(name = "Table S4 - PL a", dep = "Escape distance", model = pss, nsim = 5000)
  psx_out <- m_out(name = "Table S4 - PL b", dep = "Escape distancey", model = psx, nsim = 5000)

  out_FID_s <- rbind(mss_out, msx_out, fss_out, fsx_out, pss_out, psx_out, css_out, csx_out, hss_out, hsx_out, ass_out, asx_out, fill = TRUE)
  out_FID_s[is.na(out_FID_s)] <- ""
  out_FID_s$R2_mar = out_FID_s$R2_con = NULL
  out_FID_s[, effect := gsub("scale\\(StringencyIndex\\)", "stringency index", effect)]
  out_FID_s[, effect := gsub("scale\\(Year\\)", "year", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
  out_FID_s[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
  out_FID_s[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
  out_FID_s[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
  out_FID_s[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
  out_FID_s[, effect := gsub("Species", "species", effect)]
  out_FID_s[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_s[, effect := gsub("IDLocality", "site", effect)]
  out_FID_s[, effect := gsub("sp_loc", "species within site", effect)]
  out_FID_s[type == "random" & grepl("stringency index", effect, fixed = TRUE), effect := paste("stringency index (slope) |", gsub(" stringency index", "", effect))]
  fwrite(file = here::here("Outputs/Table_S4_rev.csv"), out_FID_s)

  load(here::here("Data/dat_est_Stringency_rev.Rdata"))
  os[predictor %in% c("scale(StringencyIndex)"), predictor := "Stringency Index"]
  oso <- os[predictor %in% c("Stringency Index")]
  oso[, N:=as.numeric(sub('.*N = ', '', model))]
  # add meta-analytical mean
    oso_s = oso[control_for_starting_distance == 'yes']
    met = summary(meta.summaries(d = oso_s$estimate, se = oso_s$sd, method = "fixed", weights = oso_s$N))$summci
    oso_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

    oso_sx = oso[control_for_starting_distance == "no"]
    metx = summary(meta.summaries(d = oso_sx$estimate, se = oso_sx$sd, method = "fixed", weights = oso_sx$N))$summci
    oso_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
    
    oso = rbind(oso, oso_met, oso_metx)
      
  oso[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

  # prepare for adding N
  oso[, N := as.character(N)]
  oso[control_for_starting_distance == "no" | is.na(N), N := ""]
  oso[, n_pos := 0.35]

  width_ <- .5 # spacing between error bars

  #col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
  #Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
  #Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
  #Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

  # From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
  #Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
  #col_ = Okabe_Ito[7:1]
  # JAMA and LocusZoom modified order
  #col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
  #col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
  col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
  #show_col(col_)
  gs6a = 
  ggplot(oso, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
      geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
      geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
      # geom_point(position = ggstance::position_dodgev(.6)) +
      geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
      # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
      # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
      # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
      # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
      geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
      scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
      #scale_color_jama(guide = "none")+ #, palette = 'light'
      scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
      scale_x_continuous(breaks = round(seq(-0.3, 0.4, by = 0.1), 1)) +
      ylab("") +
      xlab("Standardized effect size of\nStringency Index\n[on flight initiation distance]") +
      labs(tag = 'a)')+
      # coord_cartesian(xlim = c(-.15, .15)) +
      # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
      theme_bw() +
      theme(
          legend.position = "right",
          plot.tag = element_text(size = 7),
          legend.title = element_text(size = 7),
          legend.text = element_text(size = 6),
          # legend.spacing.y = unit(0.1, 'cm'),
          legend.key.height = unit(0.5, "line"),
          legend.margin = margin(0, 0, 0, 0),
          # legend.position=c(0.5,1.6),
          plot.title = element_text(color = "grey", size = 7),
          plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 1, unit = "pt"),
          panel.grid = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = ax_lines, size = 0.25),
          axis.line.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
          # axis.text.x = element_text()
          axis.ticks.length = unit(1, "pt"),
          axis.text.x = element_text(, size = 6),
          axis.text.y = element_text(colour = "black", size = 7),
          axis.title = element_text(size = 7)
      )

  if(save_plot==TRUE){
  ggsave(here::here("Outputs/Fig_S6a_Stringency.png"), gs6a, width = 8, height = 6.5, unit = "cm", dpi = 600)
  }
# predictions for Fig and Table - Google Mobility
 # full
 mgs <- lmer(scale(log(FID)) ~
     scale(Year) +
     scale(log(SD)) +
     scale(log(FlockSize)) +
     scale(log(BodyMass)) +
     scale(sin(rad)) + scale(cos(rad)) +
     # scale(Day)+
     scale(Temp) +
     scale(parks_percent_change_from_baseline) +
     (1|weekday) +  (1| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline)| Country) + (1 | IDLocality) + (1 | sp_loc),
 data = ss, REML = FALSE,
 control = lmerControl(
     optimizer = "optimx", optCtrl = list(method = "nlminb")
 )
 )
 est_mgs <- est_out(mgs, "ALL: (1|weekday) + (1|genus) + (1|Species)  + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
 est_mgs[, control_for_starting_distance := "yes"]

 mgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    #scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(StringencyIndex) +
    scale(parks_percent_change_from_baseline) +
     (1|weekday) + (1| genus) + (1 | Species) + (1 | sp_day_year) + 
     (scale(parks_percent_change_from_baseline) | Country) + (1 | IDLocality) + (1 | sp_loc),
 data = ss, REML = FALSE,
 control = lmerControl(
    optimizer = "optimx", optCtrl = list(method = "nlminb")
)
)
est_mgx <- est_out(mgx, "ALL: (1|weekday) + (1|genus) + (1|Species) + (1|sp_day_year) + (scale(Google)|Country) + (1|IDLocality) +(1|sp_loc)")
est_mgx[, control_for_starting_distance := "no"]


# CZ - singular fits only due to genera estimated as zero (removing it changes no results)
  cgs <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = ss[Country == "Czechia"], REML = FALSE
    )
  est_cgs <- est_out(cgs, "Czechia:(1|weekday)+(1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_cgs[, control_for_starting_distance := "yes"]
  
  cgx <- lmer(scale(log(FID)) ~
    #scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
    #(1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
    data = ss[Country == "Czechia"], REML = FALSE
    )
  est_cgx <- est_out(cgx, "Czechia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
  est_cgx[, control_for_starting_distance := "no"]

# FI
fgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgs <- est_out(fgs, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fgs[, control_for_starting_distance := "yes"]

fgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1| Species) + (1 | sp_day_year) + (1| IDLocality) + (1|sp_loc),
data = ss[Country == "Finland"], REML = FALSE
)
est_fgx <- est_out(fgx, "Finland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(scale(StringencyIndex)|IDLocality)+(1|sp_loc)")
est_fgx[, control_for_starting_distance := "no"]

# HU - singular fits only due to sp_loc estimated as zero (removing it changes no results)
hgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
     (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgs <- est_out(hgs, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgs[, control_for_starting_distance := "yes"]

hgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Hungary"], REML = FALSE
)
est_hgx <- est_out(hgx, "Hungary: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_hgx[, control_for_starting_distance := "no"]

# AU - singular fits only due to Year and random slope estimated as zero (removing those changes no results)
ags <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
      (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], REML = FALSE
)
est_ags <- est_out(ags, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_ags[, control_for_starting_distance := "yes"]

agx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality) + (1 | sp_loc),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Australia"], , REML = FALSE,
        control = lmerControl( 
            optimizer ='optimx', optCtrl=list(method='nlminb')) 
)
est_agx <- est_out(agx, "Australia: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)+(1|IDLocality)+(1|sp_loc)")
est_agx[, control_for_starting_distance := "no"]

# PL 
pgs <- lmer(scale(log(FID)) ~
    scale(Year) +
    scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1|genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgs <- est_out(pgs, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgs[, control_for_starting_distance := "yes"]

pgx <- lmer(scale(log(FID)) ~
    scale(Year) +
    # scale(log(SD)) +
    scale(log(FlockSize)) +
    scale(log(BodyMass)) +
    scale(sin(rad)) + scale(cos(rad)) +
    # scale(Day)+
    scale(Temp) +
    scale(parks_percent_change_from_baseline) +
    (1|weekday) + (1 | genus) + (1 | Species)+ (1 | sp_day_year),
# (1 | Year) + (1 | weekday) + (1 | genus) + (1 | Species) + (1 | sp_day_year) + (1 | IDLocality),
data = ss[Country == "Poland"], REML = FALSE
)
est_pgx <- est_out(pgx, "Poland: (1|weekday) + (1|genus)+(1|Species)+(1|sp_day_year)")
est_pgx[, control_for_starting_distance := "no"]

  # combine
    est_mgs[, Country := 'All\n(mixed model)']
    est_mgx[, Country := "All\n(mixed model)"]
    est_ags[, Country := "Australia"]
    est_agx[, Country := "Australia"]
    est_cgs[, Country := "Czechia"]
    est_cgx[, Country := "Czechia"]
    est_hgs[, Country := "Hungary"]
    est_hgx[, Country := "Hungary"]
    est_pgs[, Country := "Poland"]
    est_pgx[, Country := "Poland"]
    est_fgs[, Country := "Finland"]
    est_fgx[, Country := "Finland"]

    og = rbind(est_mgs, est_mgx, 
            est_ags, est_agx, 
            est_cgs, est_cgx, 
            est_hgs, est_hgx,
            est_pgs, est_pgx, 
            est_fgs, est_fgx)
    save(og, file = here::here('Data/dat_est_Google_rev.Rdata'))
  # estimatees for table
  mgs_out <- m_out(name = "Table S5 - full a", dep = "Escape distance", model = mgs, nsim = 5000)
  mgx_out <- m_out(name = "Table S5 - full b", dep = "Escape distance", model = mgx, nsim = 5000)
  cgs_out <- m_out(name = "Table S5 - CZ a", dep = "Escape distance", model = cgs, nsim = 5000)
  cgx_out <- m_out(name = "Table S5 - CZ b", dep = "Escape distance", model = cgx, nsim = 5000)
  fgs_out <- m_out(name = "Table S5 - FI a", dep = "Escape distance", model = fgs, nsim = 5000)
  fgx_out <- m_out(name = "Table S5 - FI b", dep = "Escape distance", model = fgx, nsim = 5000)
  hgs_out <- m_out(name = "Table S5 - HU a", dep = "Escape distance", model = hgs, nsim = 5000)
  hgx_out <- m_out(name = "Table S5 - HU b", dep = "Escape distance", model = hgx, nsim = 5000)
  ags_out <- m_out(name = "Table S5 - AU a", dep = "Escape distance", model = ags, nsim = 5000)
  agx_out <- m_out(name = "Table S5 - AU b", dep = "Escape distancey", model = agx, nsim = 5000)
  pgs_out <- m_out(name = "Table S5 - PL a", dep = "Escape distance", model = pgs, nsim = 5000)
  pgx_out <- m_out(name = "Table S5 - PL b", dep = "Escape distancey", model = pgx, nsim = 5000)

  out_FID_g <- rbind(mgs_out, mgx_out, fgs_out, fgx_out, pgs_out, pgx_out, cgs_out, cgx_out, hgs_out, hgx_out, ags_out, agx_out, fill = TRUE)
  out_FID_g[is.na(out_FID_g)] <- ""
  out_FID_g[, effect := gsub("scale\\(parks_percent_change_from_baseline\\)", "Google Mobility", effect)]
  out_FID_g[, effect := gsub("scale\\(Year\\)", "year", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(SD\\)", "starting distance (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(Temp\\)", "temperaturre", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(FlockSize\\)\\)", "flock size (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(log\\(BodyMass\\)\\)", "body mass (ln)", effect)]
  out_FID_g[, effect := gsub("scale\\(sin\\(rad\\)\\)", "time (sine of radians)", effect)]
  out_FID_g[, effect := gsub("scale\\(cos\\(rad\\)\\)", "time (cosine of radians)", effect)]
  out_FID_g[, effect := gsub("Species", "species", effect)]
  out_FID_g[, effect := gsub("sp_day_year", "species within day & year", effect)]
  out_FID_g[, effect := gsub("IDLocality", "site", effect)]
  out_FID_g[, effect := gsub("sp_loc", "species within site", effect)]
  out_FID_s[type == "random" & grepl("Google Mobility", effect, fixed = TRUE), effect := paste("Google Mobility (slope) |", gsub(" Google Mobility", "", effect))]
  fwrite(file = here::here("Outputs/Table_S5_rev.csv"), out_FID_g)

load(here::here("Data/dat_est_Google_rev.Rdata"))
og[predictor %in% c("scale(parks_percent_change_from_baseline)"), predictor := "Google Mobility"]
ogo <- og[predictor %in% c("Google Mobility")]
ogo[, N:=as.numeric(sub('.*N = ', '', model))]
# add meta-analytical mean
  ogo_s = ogo[control_for_starting_distance == 'yes']
  met = summary(meta.summaries(d = ogo_s$estimate, se = ogo_s$sd, method = "fixed", weights = ogo_s$N))$summci
  ogo_met = data.table(predictor = "Period", estimate = met[2], lwr = met[1], upr = met[3], sd = NA, model = NA, control_for_starting_distance = "yes", Country = "Combined\n(metanalytical)", N = NA)

  ogo_sx = ogo[control_for_starting_distance == "no"]
  metx = summary(meta.summaries(d = ogo_sx$estimate, se = ogo_sx$sd, method = "fixed", weights = ogo_sx$N))$summci
  ogo_metx = data.table(predictor = "Period", estimate = metx[2], lwr = metx[1], upr = metx[3], sd = NA, model = NA, control_for_starting_distance = "no", Country = "Combined\n(metanalytical)", N = NA)
  
  ogo = rbind(ogo, ogo_met, ogo_metx)
    
ogo[, Country := factor(Country, levels = rev(c("Finland", "Poland", "Czechia", "Hungary", "Australia", "Combined\n(metanalytical)", "All\n(mixed model)")))]

# prepare for adding N
ogo[, N := as.character(N)]
ogo[control_for_starting_distance == "no" | is.na(N), N := ""]
ogo[, n_pos := .15]

width_ <- .5 # spacing between error bars

#col_ <- c(brewer.pal(n = 12, name = "Paired"), "grey30", "grey80")
#Tol_bright <- c("#EE6677", "#228833", "#4477AA", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
#Tol_muted <- c("#88CCEE", "#44AA99", "#117733", "#332288", "#DDCC77", "#999933", "#CC6677", "#882255", "#AA4499", "#DDDDDD")
#Tol_light <- c("#BBCC33", "#AAAA00", "#77AADD", "#EE8866", "#EEDD88", "#FFAABB", "#99DDFF", "#44BB99", "#DDDDDD")

# From Color Universal Design (CUD): https://jfly.uni-koeln.de/color/
#Okabe_Ito <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000")
#col_ = Okabe_Ito[7:1]
# JAMA and LocusZoom modified order
#col_ =  c("#374E55FF", "#374E55FF", "#DF8F44FF", "#79AF97FF", "#00A1D5FF", "#B24745FF",  "#80796BFF") #"#6A6599FF",
#col_ <- c("#357EBDFF", "#9632B8FF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#D43F3AFF", "#D43F3AFF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
col_ = c("#357EBDFF", "#D43F3AFF", "#46B8DAFF", "#5CB85CFF", "#EEA236FF", "#9632B8FF", "#9632B8FF")[7:1] # "#D43F3AFF", "#B8B8B8FF"
#show_col(col_)

gs6b = 
ggplot(ogo, aes(x = estimate, y = Country, col = Country, shape = control_for_starting_distance)) +
    geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
    geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
    # geom_point(position = ggstance::position_dodgev(.6)) +
    geom_point(position = position_dodge(width = width_), bg = "white", size = 1.1) +
    # scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
    # scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
    # geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
    # annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
    geom_text( aes(x = n_pos,label = N), vjust = 1, size = 1.75, position = ggstance::position_dodgev(width_))+
    scale_shape_manual(name = "Controlled for\nstarting distance", guide = guide_legend(reverse = TRUE), values = c(21, 19)) +
    #scale_color_jama(guide = "none")+ #, palette = 'light'
    scale_color_manual(guide = "none", values = col_) + #guide_legend(reverse = TRUE)
    scale_x_continuous(breaks = round(seq(-0.3, 0.2, by = 0.1), 1)) +
    ylab("") +
    xlab("Standardized effect size of\nGoogle Mobility (human presence)\n[on flight initiation distance]") +
    labs(tag = 'b)')+
    # coord_cartesian(xlim = c(-.15, .15)) +
    # scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
    theme_bw() +
    theme(
        plot.tag = element_text(size = 7),
        legend.position = "right",
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 6),
        # legend.spacing.y = unit(0.1, 'cm'),
        legend.key.height = unit(0.5, "line"),
        legend.margin = margin(0, 0, 0, 0),
        # legend.position=c(0.5,1.6),
        plot.title = element_text(color = "grey", size = 7),
        plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 1, unit = "pt"),
        panel.grid = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = ax_lines, size = 0.25),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
        # axis.text.x = element_text()
        axis.ticks.length = unit(1, "pt"),
        axis.text.x = element_text(, size = 6),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 7)
    )
if(save_plot==TRUE){
ggsave(here::here("Outputs/Fig_S6b_Google_rev_width_CustomLocusZoom_v2.png"), gs6b, width = 8, height = 6.5, unit = "cm", dpi = 600)
}
#gg_S6 <- ggarrange(
 # gs6a+theme(legend.position = "none"), gs6b,
  #nrow = 1, widths=c(1, 1.3)
#)
gg_S6 <- ggarrange(
  gs6a, NULL, gs6b, widths = c(1,0.05,1),
  nrow = 1, common.legend = TRUE, legend = 'right'
)

if (save_plot == TRUE) {
  ggsave(here::here("Outputs/Fig_S6_rev.png"), gg_S6, width = 8*1.5, height = 6.5, unit = "cm", dpi = 600)
}
gg_S6

Figure S6 | Changes in avian tolerance towards humans in response to (a) stringency of governmental measures, and (b) Google Mobility. The dots with horizontal lines represent estimated standardised effect size and their 95% confidence intervals, the numbers sample sizes. For the country-specific and „All“, the effect sizes and 95% confidence intervals come from the joint posterior distribution of 5000 simulated values generated by the sim function from the arm package (Gelman et al., 2022) using the mixed model outputs controlled for starting distance of the observer (filled circles) or not (empty circles; Table S4 and S5). The models were further controlled for year, flock size (ln-transfomred), body size (ln-transformed), temperature (also a proxy for a day within the breeding season: r Pearson = 0.48; Fig. S2), and time of a day, as well as for the non-independence of data points by fitting random intercepts of weekday, genus, species, species at a given day and year, country (in All - a global mixed model), site, and species within a site, while fitting, in case of ‘All’ (global model an all countries) stringency index or Google Mobility as a random slope within Country (i.e. allowing for different effect at each Country). Fitting stringency index or Google Mobility as random slope at other random intercepts produces similar results (Fig. S1b and S1c, Table S2b and S2c). The multicollinearity was small as correlations between predictors were weak (Fig. S2). For the “Combined“, the estimate and 95% confidence interval represent the meta-analytical mean based on the country estimates and their standard deviation (from the country-specific models), and sample size per country. These analyses were restricted to data collected in the period during the COVID-19 shutdowns. Stringency index data were sourced from (Hale et al. 2021), Google Mobility from Google Mobility Reports. Depicted are estimates from country- specific models, meta-analytical mean (estimated using the country estimates, their standard deviation, and sample size per country), and a full mixed model with all countries included. The full model was controlled for starting distance (ln-transformed; filled circles) or not (empty circles), year, flock size (ln-transformed), body size (ln-transformed), temperature, and daytime. We accounted for the non-independence of data points by fitting random intercepts of weekday, genus, species, species at a given day and year, country, site, and species within a site. All continuous variables were standardised by subtracting the mean and dividing by the standard deviation.

Table S4 | Escape distance in relations to stringency index

out_FID_s$error_structure = out_FID_s$response = NULL
out_FID_s[model != "", model := c(
  "All countries", "All countries, without starting distance",
  "Finland", "Finland, without starting distance",
  "Poland", "Poland, without starting distance",
  "Czechia", "Czechia, without starting distance",
  "Hungary", "Hungary, without starting distance",
  "Australia", "Australia, without starting distance"
)]
setnames(out_FID_s, old = c("estimate_r", "lwr_r", "upr_r"), new = c("estimate", "lower", "upper"))
out_FID_s %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model N type effect estimate lower upper
All countries 3676 fixed (Intercept) -0.08 -0.332 0.176
fixed year 0.019 -0.028 0.068
fixed starting distance (ln)) 0.491 0.46 0.523
fixed flock size (ln) -0.003 -0.029 0.024
fixed body mass (ln) 0.022 -0.046 0.091
fixed time (sine of radians) -0.021 -0.062 0.021
fixed time (cosine of radians) -0.005 -0.04 0.032
fixed temperaturre -0.007 -0.05 0.034
fixed stringency index 0.017 -0.15 0.182
random species within day & year (Intercept) 5% 4% 7%
random species within site (Intercept) 5% 4% 6%
random species (Intercept) 8% 6% 8%
random site (Intercept) 8% 7% 9%
random genus (Intercept) 4% 4% 4%
random Country (Intercept) 1% -14% 14%
random stringency index (slope) | Country 1% -14% 14%
random Residual 68% 48% 93%
All countries, without starting distance 3676 fixed (Intercept) -0.113 -0.449 0.223
fixed year 0.086 0.031 0.139
fixed flock size (ln) 0.01 -0.02 0.04
fixed body mass (ln) 0.135 0.048 0.223
fixed time (sine of radians) 0.045 -0.003 0.093
fixed time (cosine of radians) 0.049 0.008 0.088
fixed temperaturre 0.009 -0.038 0.056
fixed stringency index 0.011 -0.217 0.246
random species within day & year (Intercept) 6% 4% 10%
random species within site (Intercept) 6% 4% 9%
random species (Intercept) 10% 7% 14%
random site (Intercept) 11% 8% 15%
random genus (Intercept) 6% 5% 8%
random Country (Intercept) 1% -27% 18%
random stringency index (slope) | Country 1% -27% 18%
random Residual 59% 36% 98%
Finland 354 fixed (Intercept) 0.03 -0.176 0.233
fixed year -0.108 -0.346 0.119
fixed starting distance (ln)) 0.188 0.087 0.292
fixed flock size (ln) -0.111 -0.208 -0.009
fixed body mass (ln) 0.175 0.013 0.338
fixed time (sine of radians) 0.06 -0.065 0.188
fixed time (cosine of radians) 0.012 -0.114 0.132
fixed temperaturre -0.084 -0.196 0.029
fixed stringency index 0.153 -0.066 0.384
random species within site (Intercept) 12% 12% 13%
random species within day & year (Intercept) 3% 3% 3%
random site (Intercept) 12% 10% 15%
random species (Intercept) 0% 0% 1%
random genus (Intercept) 3% 2% 5%
random Residual 68% 63% 73%
Finland, without starting distance 354 fixed (Intercept) -0.011 -0.224 0.212
fixed year -0.066 -0.298 0.167
fixed flock size (ln) -0.103 -0.202 -0.004
fixed body mass (ln) 0.245 0.082 0.409
fixed time (sine of radians) 0.056 -0.073 0.187
fixed time (cosine of radians) 0.007 -0.121 0.129
fixed temperaturre -0.085 -0.2 0.029
fixed stringency index 0.139 -0.088 0.372
random species within site (Intercept) 14% 14% 14%
random species within day & year (Intercept) 4% 3% 4%
random site (Intercept) 14% 11% 16%
random species (Intercept) 3% 2% 4%
random genus (Intercept) 1% 1% 2%
random Residual 64% 59% 69%
Poland 329 fixed (Intercept) 0.004 -0.13 0.135
fixed year -0.159 -0.391 0.068
fixed starting distance (ln)) 0.533 0.446 0.62
fixed flock size (ln) -0.024 -0.104 0.059
fixed body mass (ln) -0.095 -0.204 0.016
fixed time (sine of radians) 0.104 -0.153 0.359
fixed time (cosine of radians) -0.156 -0.419 0.102
fixed temperaturre -0.172 -0.277 -0.066
fixed stringency index 0.027 -0.196 0.255
random species within day & year (Intercept) 13% 12% 14%
random species (Intercept) 8% 5% 10%
random genus (Intercept) 0% 0% 0%
random Residual 79% 76% 82%
Poland, without starting distance 329 fixed (Intercept) 0.001 -0.201 0.207
fixed year -0.363 -0.631 -0.097
fixed flock size (ln) 0.062 -0.035 0.159
fixed body mass (ln) -0.028 -0.182 0.125
fixed time (sine of radians) -0.008 -0.311 0.31
fixed time (cosine of radians) -0.037 -0.364 0.264
fixed temperaturre -0.162 -0.287 -0.042
fixed stringency index -0.028 -0.304 0.238
random species within day & year (Intercept) 7% 7% 7%
random species (Intercept) 22% 17% 26%
random genus (Intercept) 0% 0% 0%
random Residual 72% 67% 77%
Czechia 1054 fixed (Intercept) 0.226 0.019 0.431
fixed starting distance (ln)) 0.354 0.293 0.415
fixed flock size (ln) 0.005 -0.047 0.054
fixed body mass (ln) 0.183 0.025 0.336
fixed time (sine of radians) 0.062 -0.013 0.138
fixed time (cosine of radians) 0.006 -0.065 0.076
fixed temperaturre 0.009 -0.078 0.098
fixed stringency index -0.062 -0.183 0.055
random species within site (Intercept) 4% 4% 4%
random species within day & year (Intercept) 2% 2% 2%
random species (Intercept) 15% 12% 17%
random genus (Intercept) 14% 11% 18%
random site (Intercept) 5% 3% 6%
random Residual 61% 54% 68%
Czechia, without starting distance 1054 fixed (Intercept) 0.226 0.017 0.432
fixed starting distance (ln)) 0.354 0.297 0.414
fixed flock size (ln) 0.003 -0.047 0.054
fixed body mass (ln) 0.184 0.024 0.344
fixed time (sine of radians) 0.062 -0.012 0.137
fixed time (cosine of radians) 0.006 -0.064 0.078
fixed temperaturre 0.009 -0.078 0.098
fixed stringency index -0.065 -0.182 0.054
random species within site (Intercept) 4% 4% 4%
random species within day & year (Intercept) 2% 2% 2%
random species (Intercept) 15% 12% 17%
random genus (Intercept) 14% 10% 18%
random site (Intercept) 5% 3% 6%
random Residual 61% 54% 68%
Hungary 874 fixed (Intercept) 0.073 -0.193 0.355
fixed year 0.194 0.083 0.302
fixed starting distance (ln)) 0.508 0.446 0.57
fixed flock size (ln) 0.002 -0.049 0.056
fixed body mass (ln) 0 -0.14 0.142
fixed time (sine of radians) -0.123 -0.191 -0.057
fixed time (cosine of radians) 0.034 -0.037 0.102
fixed temperaturre -0.009 -0.101 0.088
fixed stringency index 0.045 -0.092 0.187
random species within day & year (Intercept) 3% 3% 3%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 1% 1% 1%
random genus (Intercept) 11% 8% 15%
random site (Intercept) 10% 7% 13%
random Residual 74% 67% 81%
Hungary, without starting distance 874 fixed (Intercept) 0.071 -0.407 0.533
fixed year 0.24 0.112 0.368
fixed flock size (ln) 0.004 -0.053 0.063
fixed body mass (ln) 0.293 0.073 0.52
fixed time (sine of radians) -0.14 -0.219 -0.063
fixed time (cosine of radians) 0.032 -0.048 0.115
fixed temperaturre -0.031 -0.146 0.087
fixed stringency index 0.006 -0.153 0.167
random species within day & year (Intercept) 5% 5% 6%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 1% 1% 1%
random genus (Intercept) 25% 17% 31%
random site (Intercept) 18% 16% 19%
random Residual 51% 44% 60%
Australia 1065 fixed (Intercept) -0.048 -0.186 0.09
fixed year 0.036 -0.029 0.1
fixed starting distance (ln)) 0.495 0.445 0.548
fixed flock size (ln) 0.024 -0.024 0.07
fixed body mass (ln) -0.054 -0.149 0.042
fixed time (sine of radians) -0.036 -0.108 0.035
fixed time (cosine of radians) -0.052 -0.119 0.017
fixed temperaturre 0.032 -0.028 0.096
fixed stringency index 0.059 -0.004 0.122
random species within day & year (Intercept) 6% 6% 6%
random species within site (Intercept) 12% 12% 12%
random species (Intercept) 11% 9% 13%
random genus (Intercept) 2% 2% 3%
random site (Intercept) 7% 5% 8%
random Residual 61% 57% 66%
Australia, without starting distance 1065 fixed (Intercept) -0.101 -0.279 0.083
fixed year 0.101 0.027 0.176
fixed flock size (ln) 0.054 0 0.11
fixed body mass (ln) 0.06 -0.065 0.181
fixed time (sine of radians) -0.01 -0.095 0.074
fixed time (cosine of radians) -0.021 -0.096 0.055
fixed temperaturre 0.039 -0.033 0.11
fixed stringency index 0.09 0.016 0.163
random species within day & year (Intercept) 7% 7% 8%
random species within site (Intercept) 11% 11% 11%
random species (Intercept) 13% 11% 15%
random genus (Intercept) 7% 6% 8%
random site (Intercept) 8% 6% 10%
random Residual 54% 49% 58%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance.

Table S5 | Escape distance in relations to Google Mobilitty

out_FID_g$R2_mar = out_FID_g$R2_con = out_FID_g$error_structure = out_FID_g$response = NULL
out_FID_g[model != "", model := c(
  "All countries", "All countries, without starting distance",
  "Finland", "Finland, without starting distance",
  "Poland", "Poland, without starting distance",
  "Czechia", "Czechia, without starting distance",
  "Hungary", "Hungary, without starting distance",
  "Australia", "Australia, without starting distance"
)]
setnames(out_FID_g, old = c('estimate_r','lwr_r','upr_r'), new =c('estimate','lower','upper' ))
out_FID_g %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
model N type effect estimate lower upper
All countries 3644 fixed (Intercept) 0.02 -0.096 0.136
fixed year 0.016 -0.024 0.059
fixed starting distance (ln)) 0.499 0.466 0.53
fixed flock size (ln) -0.001 -0.027 0.025
fixed body mass (ln) 0.026 -0.039 0.092
fixed time (sine of radians) -0.009 -0.049 0.03
fixed time (cosine of radians) 0.011 -0.022 0.044
fixed temperaturre 0.022 -0.015 0.061
fixed Google Mobility -0.096 -0.225 0.034
random species within day & year (Intercept) 6% 5% 7%
random species within site (Intercept) 5% 4% 5%
random species (Intercept) 7% 7% 7%
random site (Intercept) 9% 9% 9%
random genus (Intercept) 0% -1% 4%
random genus Google Mobility 0% -1% 4%
random Country (Intercept) 1% -4% 4%
random Country Google Mobility 1% -4% 4%
random Residual 71% 59% 82%
All countries, without starting distance 3644 fixed (Intercept) 0.111 -0.133 0.361
fixed year 0.079 0.029 0.129
fixed flock size (ln) 0.012 -0.017 0.042
fixed body mass (ln) 0.133 0.047 0.218
fixed time (sine of radians) 0.067 0.02 0.115
fixed time (cosine of radians) 0.075 0.037 0.113
fixed temperaturre 0.053 0.004 0.1
fixed scale(StringencyIndex) 0.047 -0.026 0.118
fixed Google Mobility -0.116 -0.296 0.066
random species within day & year (Intercept) 7% 5% 9%
random species within site (Intercept) 6% 4% 8%
random species (Intercept) 10% 8% 12%
random site (Intercept) 13% 10% 15%
random genus (Intercept) 0% -1% 5%
random genus Google Mobility 0% -1% 5%
random Country (Intercept) 1% -16% 11%
random Country Google Mobility 1% -16% 11%
random Residual 62% 41% 88%
Finland 322 fixed (Intercept) 0.022 -0.187 0.223
fixed year 0.132 -0.007 0.268
fixed starting distance (ln)) 0.194 0.084 0.295
fixed flock size (ln) -0.119 -0.222 -0.013
fixed body mass (ln) 0.243 0.084 0.398
fixed time (sine of radians) 0.132 -0.009 0.269
fixed time (cosine of radians) 0.041 -0.086 0.167
fixed temperaturre -0.004 -0.139 0.13
fixed Google Mobility -0.102 -0.241 0.041
random species within site (Intercept) 13% 12% 13%
random species within day & year (Intercept) 3% 2% 3%
random site (Intercept) 17% 13% 20%
random species (Intercept) 0% 0% 0%
random genus (Intercept) 2% 1% 3%
random Residual 66% 61% 71%
Finland, without starting distance 322 fixed (Intercept) -0.008 -0.227 0.213
fixed year 0.152 0.01 0.295
fixed flock size (ln) -0.111 -0.218 -0.007
fixed body mass (ln) 0.321 0.163 0.479
fixed time (sine of radians) 0.12 -0.031 0.266
fixed time (cosine of radians) 0.031 -0.094 0.158
fixed temperaturre -0.02 -0.153 0.116
fixed Google Mobility -0.088 -0.236 0.057
random species within site (Intercept) 14% 14% 15%
random species within day & year (Intercept) 3% 3% 3%
random site (Intercept) 18% 14% 21%
random species (Intercept) 0% 0% 0%
random genus (Intercept) 2% 1% 4%
random Residual 62% 57% 67%
Poland 329 fixed (Intercept) 0.004 -0.13 0.137
fixed year -0.19 -0.324 -0.049
fixed starting distance (ln)) 0.533 0.447 0.616
fixed flock size (ln) -0.023 -0.105 0.06
fixed body mass (ln) -0.094 -0.204 0.011
fixed time (sine of radians) 0.109 -0.13 0.358
fixed time (cosine of radians) -0.163 -0.418 0.083
fixed temperaturre -0.174 -0.273 -0.079
fixed Google Mobility 0.007 -0.119 0.128
random species within day & year (Intercept) 13% 12% 14%
random species (Intercept) 7% 5% 10%
random genus (Intercept) 0% 0% 0%
random Residual 79% 76% 82%
Poland, without starting distance 329 fixed (Intercept) 0.004 -0.202 0.211
fixed year -0.379 -0.532 -0.225
fixed flock size (ln) 0.061 -0.041 0.158
fixed body mass (ln) -0.027 -0.178 0.129
fixed time (sine of radians) -0.011 -0.296 0.275
fixed time (cosine of radians) -0.041 -0.325 0.259
fixed temperaturre -0.162 -0.28 -0.045
fixed Google Mobility 0.053 -0.086 0.193
random species within day & year (Intercept) 7% 7% 7%
random species (Intercept) 21% 16% 26%
random genus (Intercept) 0% 0% 0%
random Residual 72% 67% 77%
Czechia 1054 fixed (Intercept) 0.229 0.018 0.431
fixed starting distance (ln)) 0.345 0.286 0.403
fixed flock size (ln) 0.009 -0.041 0.059
fixed body mass (ln) 0.186 0.026 0.345
fixed time (sine of radians) 0.063 -0.006 0.134
fixed time (cosine of radians) 0.034 -0.019 0.089
fixed temperaturre 0.066 -0.006 0.136
fixed Google Mobility -0.09 -0.161 -0.018
random species within site (Intercept) 5% 4% 5%
random species within day & year (Intercept) 1% 0% 1%
random species (Intercept) 14% 12% 17%
random genus (Intercept) 14% 10% 18%
random site (Intercept) 5% 4% 6%
random Residual 61% 54% 69%
Czechia, without starting distance 1054 fixed (Intercept) 0.226 0.022 0.432
fixed starting distance (ln)) 0.344 0.284 0.402
fixed flock size (ln) 0.009 -0.04 0.059
fixed body mass (ln) 0.184 0.026 0.342
fixed time (sine of radians) 0.063 -0.007 0.132
fixed time (cosine of radians) 0.035 -0.024 0.09
fixed temperaturre 0.066 -0.007 0.139
fixed Google Mobility -0.09 -0.163 -0.016
random species within site (Intercept) 5% 4% 5%
random species within day & year (Intercept) 1% 0% 1%
random species (Intercept) 14% 11% 17%
random genus (Intercept) 14% 10% 18%
random site (Intercept) 5% 4% 6%
random Residual 61% 54% 69%
Hungary 874 fixed (Intercept) 0.082 -0.182 0.351
fixed year 0.19 0.078 0.299
fixed starting distance (ln)) 0.505 0.444 0.566
fixed flock size (ln) 0.001 -0.052 0.054
fixed body mass (ln) 0.001 -0.142 0.14
fixed time (sine of radians) -0.122 -0.189 -0.055
fixed time (cosine of radians) 0.023 -0.042 0.089
fixed temperaturre -0.027 -0.106 0.051
fixed Google Mobility -0.029 -0.127 0.072
random species within day & year (Intercept) 3% 3% 3%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 1% 1% 1%
random genus (Intercept) 11% 7% 15%
random site (Intercept) 9% 6% 13%
random Residual 75% 67% 82%
Hungary, without starting distance 874 fixed (Intercept) 0.087 -0.346 0.536
fixed year 0.323 0.192 0.457
fixed flock size (ln) -0.001 -0.058 0.059
fixed body mass (ln) 0.282 0.061 0.508
fixed time (sine of radians) -0.156 -0.235 -0.076
fixed time (cosine of radians) 0.032 -0.039 0.105
fixed temperaturre -0.044 -0.136 0.05
fixed Google Mobility -0.107 -0.227 0.017
random species within day & year (Intercept) 5% 5% 5%
random species within site (Intercept) 0% 0% 0%
random species (Intercept) 1% 1% 2%
random genus (Intercept) 24% 18% 30%
random site (Intercept) 16% 13% 18%
random Residual 54% 45% 63%
Australia 1065 fixed (Intercept) -0.055 -0.194 0.091
fixed year 0.028 -0.037 0.091
fixed starting distance (ln)) 0.5 0.448 0.55
fixed flock size (ln) 0.024 -0.024 0.07
fixed body mass (ln) -0.049 -0.142 0.043
fixed time (sine of radians) -0.044 -0.12 0.03
fixed time (cosine of radians) -0.044 -0.113 0.024
fixed temperaturre 0.018 -0.039 0.077
fixed Google Mobility -0.035 -0.09 0.019
random species within day & year (Intercept) 6% 6% 6%
random species within site (Intercept) 12% 12% 12%
random species (Intercept) 12% 10% 14%
random genus (Intercept) 2% 1% 2%
random site (Intercept) 8% 6% 10%
random Residual 61% 56% 65%
Australia, without starting distance 1065 fixed (Intercept) -0.107 -0.294 0.084
fixed year 0.083 0.012 0.155
fixed flock size (ln) 0.054 -0.001 0.109
fixed body mass (ln) 0.067 -0.056 0.187
fixed time (sine of radians) -0.031 -0.114 0.056
fixed time (cosine of radians) -0.017 -0.098 0.064
fixed temperaturre 0.008 -0.061 0.075
fixed Google Mobility -0.018 -0.082 0.045
random species within day & year (Intercept) 7% 7% 7%
random species within site (Intercept) 10% 10% 11%
random species (Intercept) 14% 12% 16%
random genus (Intercept) 6% 5% 7%
random site (Intercept) 10% 8% 12%
random Residual 52% 48% 57%

Continuous variables were z-transformed. For descriptions of variables see Methods of the paper. Residual indicates residual variance.


References

  • Bulla, M., Blumstein, D.T., Benedetti, Y., Floigl, K., Jokimäki, J., Kaisanlahti-Jokimäki, M.-L., Markó, G., Morelli, F., Siretckaia, A., Szakony, S., Weston, M.A., Zeid, F.A., Tryjanowski, P., Albrecht, T. & Mikula, P. (2022). Supporting information for ‘Urban birds’ flight responses were unaffected by the COVID-19 shutdowns’. Open Science Framework https://doi.org/10.17605/OSF.IO/WUZH7.
  • Gelman, A., Su, Y.-S., Yajima, M., Hill, J., Pittau, M., Kerman, J., Zheng, T., & Vincent, D. (2016). Data Analysis using Regression and Multilevel/Hierarchical Models. In CRAN Repository (1.8-6.; pp. 1–53).
  • Hale T., Angrist N., Goldszmidt R., Kira B., Petherick A., Phillips T., Webster S., Cameron-Blake E., Hallas L., Majumdar S., Tatlow H. (2021). A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker). Nature Human Behaviour 2021 5:4 5:529–538. https://doi.org/10.1038/s41562-021-01079-8.
  • Peterson, B. G., & Carl, P. (2020). PerformanceAnalytics: Econometric Tools for Performance and Risk Analysis. R package version 2.0.4. https://CRAN.R-project.org/package=PerformanceAnalytics
  • Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag, New York.

Session info

  sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.6.2              viridisLite_0.4.1          scales_1.2.1               rphylopic_0.3.0            rmeta_3.0                 
##  [6] RColorBrewer_1.1-3         png_0.1-7                  PerformanceAnalytics_2.0.4 xts_0.12.1                 zoo_1.8-10                
## [11] performance_0.9.0          optimx_2022-4.30           multcomp_1.4-19            TH.data_1.1-1              survival_3.3-1            
## [16] mvtnorm_1.1-3              kableExtra_1.3.4           here_1.0.1                 gtable_0.3.1               ggtext_0.1.2              
## [21] ggsci_2.9                  ggpubr_0.4.0               ggimage_0.3.1              ggplot2_3.3.6              foreach_1.5.2             
## [26] effects_4.2-1              carData_3.0-5              data.table_1.14.6          arm_1.12-2                 lme4_1.1-29               
## [31] Matrix_1.4-1               MASS_7.3-57                rmarkdown_2.18            
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4         colorspace_2.0-3    ggsignif_0.6.3      rprojroot_2.0.3     gridtext_0.1.4      rstudioapi_0.14     httpcode_0.3.0     
##  [8] farver_2.1.1        fansi_1.0.3         xml2_1.3.3          codetools_0.2-18    splines_4.2.0       cachem_1.0.6        knitr_1.41         
## [15] jsonlite_1.8.4      nloptr_2.0.3        broom_0.8.0         gridBase_0.4-7      compiler_4.2.0      httr_1.4.4          backports_1.4.1    
## [22] assertthat_0.2.1    fastmap_1.1.0       survey_4.1-1        cli_3.4.1           htmltools_0.5.4     tools_4.2.0         coda_0.19-4        
## [29] glue_1.6.2          dplyr_1.0.10        Rcpp_1.0.9          jquerylib_0.1.4     vctrs_0.5.1         crul_1.2.0          svglite_2.1.0      
## [36] nlme_3.1-157        iterators_1.0.14    insight_0.18.8      xfun_0.35           stringr_1.5.0       rvest_1.0.3         lifecycle_1.0.3    
## [43] rstatix_0.7.0       sandwich_3.0-1      yaml_2.3.6          curl_4.3.3          gridExtra_2.3       ggfun_0.0.6         yulab.utils_0.0.4  
## [50] sass_0.4.4          stringi_1.7.8       highr_0.9           boot_1.3-28         rlang_1.0.6         pkgconfig_2.0.3     systemfonts_1.0.4  
## [57] evaluate_0.18       lattice_0.20-45     purrr_0.3.4         labeling_0.4.2      tidyselect_1.1.2    magrittr_2.0.3      R6_2.5.1           
## [64] magick_2.7.3        generics_0.1.3      DBI_1.1.3           pillar_1.8.1        withr_2.5.0         abind_1.4-5         nnet_7.3-17        
## [71] tibble_3.1.8        car_3.0-13          utf8_1.2.2          grid_4.2.0          digest_0.6.30       webshot_0.5.4       tidyr_1.2.0        
## [78] numDeriv_2016.8-1.1 gridGraphics_0.5-1  munsell_0.5.0       ggplotify_0.1.0     bslib_0.4.1         mitools_2.4         quadprog_1.5-8
# END